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Abstract 


An efficient and robust computational scheme is given 
for the calculation of the frequency response functions 
of a large order, flexible system implemented with a lin- 
ear, time invariant control system. Advantage is taken 
of the highly structured sparsity of the system matrices 
based on a model of the structure using normal mode co- 
ordinates. The computational time per frequency point 
of the new computational scheme is a linear function of 
system size, a significant improvement over traditional, 
full-matrix techniques whose computational times per fre- 
quency point range from quadratic to cubic functions of 
system size. This permits the practical frequency domain 
analysis of systems of much larger order than by tradi- 
tional, full-matrix techniques. Formulations are given for 
both open- and closed-loop systems. Numerical examples 
are presented showing the advantages of the present for- 
mulation over traditional approaches, both in speed and 
in accuracy. Using a model with 703 structural modes, 
the present method was up to two orders of magnitude 
faster than a traditional method. The present method gen- 
erally showed good to excellent accuracy throughout the 
range of test frequencies, while traditional methods gave 
adequate accuracy for lower frequencies, but generally 
deteriorated dramatically in performance at higher fre- 
quencies. 
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1. Introduction 


Control of flexible systems has received significant attention in the literature. 
To date, numerous techniques, algorithms and procedures have been developed for 
design of controllers for such systems ranging from spacecraft and satellites to 
aircraft, ships, machines, etc. These flexible systems which are generally infinite- 
dimensional are typically modeled using a finite number of generalized coordinates 
or modes. Control of flexible systems may become difficult depending on the 
number, location, relative proximity, and inherent damping of these modes. The 
response of the system to a given disturbance/excitation generally depends on 
modal properties (amplitude, frequency, and damping) and the amplitude and phase 
content of the disturbance/excitation. In general, two techniques, time domain 
analysis and frequency domain analysis, have been developed and extensively 
used to analyze and characterize the input/output behavior of linear time-invariant 
systems including flexible systems. In frequency domain analysis, frequency 
response functions (defined as transfer function matrices from inputs to outputs of 
the system) have typically been used (usually in the form of magnitude and phase 
or Bode plots) in the analysis of linear systems as well as in designing controllers 
for such systems. In general, frequency response functions of the open-loop system 
are used to evaluate the performance of the open-loop system, and to identify and 
quantify needed performance and/or stability improvements at various frequency 
bands. The closed-loop frequency response functions are typically needed to insure 
that desired performance and stability have been achieved by the control system. 
Moreover, frequency-domain specifications such as peak magnitude, bandwidth, 
roll-off rate, etc. are often used in characterizing the desired behavior of the system 
in the frequency domain (this is known as loop shaping). 

In general, the order of the flexible system (as defined by the number of modes 
retained in the model) for which open-loop and/or closed-loop analysis is performed 
depends on the application considered. For example, if the closed-loop response 
of a spacecraft with a low-bandwidth attitude control system is of interest, then a 
small set of modes would be sufficient to capture the low frequency closed-loop 
behavior of the system. On the other hand, if the response of the flexible system 
is desired over a large frequency range or if the control system considered has a 
high bandwidth, then a large set of modes (in the hundreds or thousands) may be 
necessary to capture the true response of the system. 
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However, the current techniques for obtaining frequency response functions, 
although able to deal with small or medium size systems, have problems in handling 
large order systems. A straightforward calculation of the frequency response 
function matrix at a single frequency point which is based on the definition of the 
transfer function has a computational cost which is a cubic function of the system 
size. If this calculation must be repeated for many frequency points, Laub ([1, 2]) 
presents a technique which has a better average cost. This technique performs an 
initial orthogonal transformation of the system which reduces the system response 
matrix to upper Hessenberg form. This initial transformation has a computational 
cost which is a cubic function of the system size. This technique can then calculate 
the frequency response function matrix at each frequency point at a cost which is 
a quadratic function of the system size. However, for very large systems (many 
hundreds of modes or more), even this is too slow, and a better method is needed. 

To this end, this paper describes efficient techniques for the computation of 
open- and closed-loop frequency response functions of large order flexible systems. 
The proposed techniques are computationally robust and accurate. The closed- 
loop frequency response function calculation is novel and constitutes enabling 
technology for the frequency domain analysis of large flexible systems. These 
techniques take advantage of the sparsity of the flexible systems in normal mode 
coordinates and reduce the computational cost from a quadratic function of the order 
of the system to a linear function. A fringe benefit of these faster computational 
techniques is an improvement in accuracy. The decoupling of the calculations 
involving each structural mode not only makes the calculation faster, but also 
means that lower order, better conditioned linear systems are being solved; and this 
improves accuracy. Numerical examples are presented showing the advantages of 
the present formulation over traditional approaches, both in speed and in accuracy. 
The present paper expands on the authors’ previous work [3]. 

Symbols 


A generic matrix 

A, k x A control system state matrix 

A s 2 p x 2 p plant state matrix for first-order model, see Equation (3) 

^ 2x2 diagonal block i of A 6 1 < / < p\ see Equation (4) 
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4 (2 p 4- r) x (2 p -f r) system matrix of the system which combines 

both the plant and the controller 

D c k x / matrix of tracking error (resp., reference input) influence on 
closed-loop (resp., open-loop) control system states 

D s 2 p x m matrix of control input influence on first-order model states, 

see Equation (5) 

B w 2 p x g matrix of disturbance input influence on first-order model 

states 

B any of the block matrix columns in Equation (31) 


b generic scalar or vector 

j(i) row i of matrix 

C c rn x k matrix of control system state influence on negative of control 

input vector 

Cp f x n matrix of influence of states of second-order physics based 

model on measurement output 

c„ f x p matrix of influence of modal model states on measurement 
output 

Cp(i,j) (/. /) element of matrix C p 

C r f x n matrix of influence of state rates of second-order physics 

based model on measurement output 

Cr f x p matrix of influence of modal model state rates on 

measurement output 

C r ( i ,j) (i.j) element of matrix C r 

C s / x 2 p matrix of influence of first-order model states on 

measurement output, see Equation (6) 

C pr I x 2 p matrix of influence of states of first-order model with input 
feedthrough on performance outputs 

C( t T / x n matrix of influence of state accelerations of second-order 
model on performance output 

C pr l x p matrix of influence of modal model state accelerations on 
performance output 
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Cf 

cV 

cV 

cr 

c pr 

r ,pr 

L 2 

c 

D 

D 

D P r 

Dfr 

D 

d 

E{*) 

En(s) 

EvM 

E 2l(«) 
E 22W 

e 

/ 

G 


l x // matrix of influence of states of second-order model on 
performance output 

/ x p matrix of influence of modal model states on performance 
output 

I x a matrix of influence of state rates of second-order model on 
performance output 

/ x p matrix of influence of modal model state rates on performance 
output 

/ x 2 P matrix of influence of states of first-order model with 
acceleration term on performance outputs 

1 x 2 j> matrix of influence of state rates of first-order model with 
acceleration term on performance outputs 
generic input influence matrix of a linear system 

n x f> damping matrix of second-order physics based model 
p x p damping matrix of modal model 

/ x m feedthrough influence of control input on performance outputs 
in first-order model with input feedthrough 
/ x y feedthrough influence of disturbance input on performance 
outputs in first-order model with input feedthrough 

generic feedthrough matrix of a linear system 
m x 1 input noise/di strubance vector 
(2 p+ k) x (2 p + k) matrix -v/ — A 
2 p x 2 p upper left submatrix of E(s) 

2 p x k upper right submatrix of E(n) 

A x 2 p lower left submatrix of E(s) 

A x A lower right submatrix of E(s) 

/ x 1 tracking error vector 

number elements in measurement output vector 

generic matrix output of Matlab function ltifr 
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G(s) (/ + /) x (>n + (j) open-loop transfer function matrix of first order 
plant model 

G'ii(.s) 1 x w upper left submatrix of G(.s-) 

G j 2 ( -s- ) I x g upper right submatrix of G'(.s) 

G' 2 i ( •*>■ ) / x in lower left submatrix of G(s) 

G-22 (••»■) / x g lower right submatrix of G(s) 

g number of elements in exogenous disturbance vector 

H it x in matrix of influence of control input in second-order model 

H p x in matrix of influence of control input in modal model 

Hjj the (i.j) element of matrix H 

H tr n x g matrix of influence of exogenous disturbance in second-order 
model 

H w v x g matrix of influence of exogenous disturbance in modal model 
I identity matrix of the proper size to make sense in context 

/ subscript to index the retained modes, 1 < / < p 

J \f— 1 , the imaginary unit for complex numbers 

A n x n stiffness matrix of second-order model 

K(a) m x / controller transfer function matrix 
K p x p stiffness matrix of modal model 

k number of control system states 

l number of elements in performance output vector 

M n x n mass matrix of second-order model 

M p x p mass matrix of modal model 

in number elements in control input vector 

n number of states in second-order model 


p number of normal modes retained for modal equations 

Q(- s ) 2 p x in matrix (si — 

Qi(s ) 2 x p submatrix of Q(s) containing rows number 21 — 1 and 2 i, 

1 < i < p 
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(j p x 1 vector of modal coordinates 

<jj component i of vector q, 1 < / < p 

r f x 1 reference input vector 

s complex variable used in transfer functions (Laplace transforms) 

T s sampling frequency of a discrete-time system 

T a j transfer function from /3 to a, where a € {/y, y pr ,p, w} and 
,3 G {r, </, u\ r} 

u m x 1 control input vector 

v / x 1 measurement noise vector 

TT' generic matrix 

w f/ x 1 exogenous disturbance vector 

X generic matrix 

A'(s) (2 ])+ k) x (2y> + A) matrix E{s)~ l 

X i j ( .s ) 2 p x 2 p upper left submatrix of X ( -s- ) 

A'i 2 («) 2px A upper right submatrix of X(v) 

A' 2 i(-v) A' x 2 j> lower left submatrix of X{s) 

X'22 ( s ) A- x A lower right submatrix of A’ (a) 

x n x 1 position/attitude vector of second-order model 

x generic scalar or vector unknown in a linear equation 

x a generic approximate scalar or vector solution to a linear equation 

x t A- x 1 vector of control system states 

x s 2 p x 1 state vector for first-order model, see Equation (2) 

Y generic matrix 

y f x 1 vector of measurements outputs 

y pr 1 x 1 vector of performance outputs 

s generic vector 

A(a) k x A- matrix £ 22 (s) - E 2 iE^ l (s)E V2 

6 symmetric relative error function 

( generic complex number 
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open-loop damping ratio of retained mode /, 1 < / < p 
generic complex number 

scalar comparison magnitude for generic scalar ./ 

n x p matrix whose columns are the p retained mode 
shapes [c>i <fo ■ ■ • ] 

n x 1 mode shape vector of retained mode i, 1 < / < p 

a scalar or vector of frequency values 

open-loop frequency of retained mode /',!</< p 




2. Mathematical Formulation 


2.1. Second-Order Modal Equations 

The dynamics of a typical linear, time-invariant flexible system may be ex- 
pressed in a second-order model as 

M:i + D.v + Kx = Hu + H w iv + Hd 
IJ = CpX + C V X 

if r = cp- + cp + cp 


where 1/, D, and K are the n x n mass, damping and stiffness matrices, respec- 
tively; x is the 77 x 1 position/attitude vector; u is the in x 1 control input vector; 
((’ is the tj x 1 exogenous disturbance vector; d is the m x 1 input noise/disturbance 
vector; H is the n x m control input influence matrix; and H w is the n x <j ex- 
ogenous disturbance influence matrix. The vectors y and y pr are, respectively, the 
/ x 1 measurements output vector and the / x 1 vector of performance outputs; c P 
and C, are / x n measurement output influence matrices; and C pr , C r 7 , and 
are / x n performance output influence matrices. 

If the second-order system is transformed into normal mode coordinates, and j> 
of the normal modes are retained to capture the relevant dynamics of the structure, 
then the system equations may be written in a modal form as 

Mq + Dq + Kq = Hu + H w w + Hd 
y = C p q + C r q 
y pT = Cp! + C? r q + Cpi 


where M, D, and K are the p x p modal mass, damping, and stiffness matrices, 
respectively; q is the p x 1 vector of modal coordinates; and H and H w are the p x m 
control input and the p x g disturbance influence matrices in modal coordinates, 
respectively. The matrices C p and C, are / x p measurement output influence 
matrices in modal coordinates; and Cf, r , C pr , and C ,pr are / x p performance 
output influence matrices in modal coordinates. 

It is assumed that the mode shapes are normalized with respect to the mass 
matrix, and modal damping is assumed. This means that il7 = /, D = 
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diag{2£iu;i, —C‘2^2 • • • • i 2 pulp ]■ , and A — • • • • 1 ^'p } where w>/ 3nd 

(; are the open-loop frequencies and damping ratios. 

The control input and disturbance influence matrices are given by: 

H = <& t H 
H w = § t H w 


The measurement and performance output influence matrices are given by: 


C p = C p * 


npr 

'-p 


; CV = C,<I> 
Cf r = cr$ ; 


CP r 
^ a 




The columns of matrix $ are the p retained mode shapes: 

$ = [on 02 ■■■ 0p] 

2.2. First-order Form of Equations 

The second-order modal equations may be rewritten in a first-order form as 

x s = A»x s -1- B s ti + B w w -I- B s d 

y = C s x s ( 1 ) 

V vr = Cfx, + Cfx,. 

The vector .r s is the plant state vector whose components are 

'<71 ' 

<71 
<72 


. Q ' 2 V 

= < • > , 


<7?> 
f <7/i 


( 2 ) 


and the vectors y and y pr are the same plant measurement and performance output 
vectors as before. The matrix is the plant state matrix and has the form 


r ,tl 


A s — 


0 

o' A\ 


0 

0 


( 3 ) 


0 0 


AS 
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where 


The matrix D # is the control input influence matrix, formed by setting its odd- 
numbered rows to zeros and using the rows of H for its even-numbered rows: 

"0 _0 O' 

H\ i H 12 H\ m 

0 _0 _0 

Hn Bn H 2 ,„ 


0 0 _0 

_Hp \ H p2 H pm _ 

in which H t] , for example, represents the ( i . j ) element of matrix H . The matrix 
B w is formed from H w in the same manner. 

The measurement output influence matrix, C 6 is defined by setting the odd- 
numbered columns of C. s to the columns of C' p and the even numbered columns 
of Co to the columns of C r : 


Cp( 1 , 1 ) 0 ( M) 

0 ( 1 , 2 ) 

0 ( 1 , 2 ) ••• 

C v (\,p) 

0 ( 1 . p) 

0 ( 2 , i) 0 ( 2 , 1 ) 

C 9 ( 2 , 2 ) 

0 ( 2 , 2 ) ... 

• * 

■■■ C p ( 2 ,p) 

0 ( 2 .p) 

C v {fA) Cr(fA) 

CpLfA) 

0 (/, 2 ) ■■■ 

Cp(Bp) 

Cr(f.p) 


Here, C p (i,j) and C r (i. j) denote the (/,./) element of matrix C p and C r , respec- 
tively. C ] " is defined from Cf/ and Cf r in the same fashion. C pr is defined 
by setting the odd numbered columns of Cf to zeros and the even numbered 
columns of Cf to the columns of C// r . 

By substituting the first equation of (1) into the third, the acceleration term can 
be replaced by feedthrough: 

i's — “I - BgU + B w w T- B s d 

V = C s . x 9 (7) 

y pr = C pr Xs + D pr u + ppr {(7 + jypr A 

The performance output influence matrix is given by 

C ,pr = Cf + 
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while the performance feedthrough matrices are 

D pr = C Pr Bs . D pr 

Notice that if there is no performance acceleration output (C'a — 0), then 
Cl” = 0 and Clf = 0, so both feedthrough matrices, D pr and D ! ”\ are zero. 


2.3. Control System Equations 

In this paper, it is assumed that the structure is controlled by a linear time- 
invariant control system. The model of a linear time-invariant control system for 
a typical flexible structure may be written as 

.r c = A c x c + D ( < 

u = C c x c (8) 

e = r — y — v 

where x c denotes the A x 1 vector of control system states; e denotes the / x 1 
tracking error vector; r denotes the / x 1 reference input vector; v is the / x 1 
measurement noise vector; A c , D c , and C\ represent the k x A control system state 
matrix, the Ax/ input influence matrix, and the rn x A output influence matrix, 
respectively; and y is the measurement output vector which was defined in the 
previous section. Note that the control system equations, as represented by Eq. 
(8), do not include any feedthrough terms. However, if there are feedthrough terms 
present in the control system, then augmented dynamics (with roll-off filters) are 
used to reduce the control system equations to the form given in Eq. (8). This 
procedure is described in [4]. If the system is open-loop, then v is zero and y is 
not fed back; so e is r, and the control system model becomes: 


x c = A c x c + Bel- 
li — C c x c 


(9) 


2.4. Equations for Open-Loop Configuration 

The block diagram of the system for the open-loop configuration is shown in 
Figure 1. 
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K(s) 


d 

u 1+ 


w 


\P T 



+ 


G(s) 


JT 


y 


Figure 1. Block diagram for open-loop configuration. 


Here, K(s) denotes the controller transfer function matrix, which is defined for 
all complex s not in the spectrum of A c as 

K(») = C c (»I-A c )-'B c (10) 


The plant transfer function matrix G(s) is partitioned according to inputs, d and u\ 
and outputs, y pr and y, and is defined for all complex s not in the spectrum of A# as 


G(h) 


G\ i (-s) G [2 (s) 

G'2 1 ( •* ) G'2'2 ( *') 


( 11 ) 


with 


G] i =CV(*/-.4 s r l B 4 + D’’ r 
G i2 = C' ,r ( 5 / - + D £, r 

G -21 = c„(.il - A S )~ : B S 
G22 = C s (sI-A s )- , B w 


( 12 ) 


The state-space models of the controller and plant, given by Eqs. (7) and (9), 
are combined to give 



•'=<> 

y pr _ y C w D pr Cc ]J^ \ + [0 D pr 

[ x c J 

„ = [° r 



(13) 

(14) 

(15) 

(16) 
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For the open-loop system model in Eqs. (13)-(16), there are three possible inputs 
to the system (i\d< ic) and three possible outputs (ij. i/ p7 \ u). Therefore, a total of 
nine transfer functions can be defined for the system, with each transfer function 
corresponding to an input/output pair. Irrespective of the choice of the input/output 
pair, the computation of open-loop transfer functions from Eqs. (13)-(16) would 
involve the inverse of matrix E{s) = (si - A^ at each frequency point * = 

where A is the system matrix of the combined plant model and the controller, 
given by 

'As BX'c 

. 0 A 

Partition E(s) conformally: 


'En(») 

E\2 


si - Ai 

-B s C c ' 

0 

El 2 (■'>■). 


0 

si - A c 


Set X(s) — E{s) 1 and partition conformally: 


r.Ynw .y I2 (s)1 


'En(s) 

Ei2 ' 

1 

yN 

'N 

r i 
>1, 


0 

E-n{*)_ 


Then, from matrix inversion lemma, one has 

X u {s) = E^{s) 

AT 12 ( ) = — E 1 1 (s) E 12 E -22 ( ) 

x 2 M = 0 

A 22 (,) = e£{s) 


(18) 


(19) 


2.5. Open-Loop Transfer Functions 

Formulae for the nine open-loop transfer functions are derived in this section. 
Each transfer function formula is stated either in terms of the plant and controller 
system matrices and terms derived from them in the previous section, or as a simple 
alteration to a previously derived transfer function. The next section will address 
issues of computational efficiency. 
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Subscripts are used in the following to relate input connections to output 
connections. The first subscript indicates the output and the second subscript 
indicates the input. 

• The transfer function from the tracking command, r, to measurement output, 
y, is given by: 


Tyr(s) = G' 2 i(.s)A» = [C, 0] 

’-1 


'A'll(v) 

A’l2(*) 

' 0 ' 

.-V 2I (») 

X 2 -2(s)_ 

p c 


( 20 ) 


The transfer function from the input noise/disturbances, d, to measurement 
output, y, is given by: 


J>W = G J1 (*) = [C, 0] 
= C.Bf/Mi?, 


A'ii(.s) A'i2(.s) 

'B,' 

.a 21 (•'>•) A 22 ( • s ' ) _ 

_ 0 _ 


( 21 ) 


The transfer function from the exogenous disturbances, w, to measurement 
output, y, is given by: 


T,A») = <? 22 <*) = C'Eu'WBr 


( 22 ) 


The transfer function from the tracking command, r, to performance output, 
is given by: 


A' u (,) 

A'l2 (*)' 

'0 ' 

|A 21 (*) 

-V 2 2(*). 

B c 


7> r (*) = Gii(*)A'(») = [C' )r D^Cc] 

= + Dr T C c X- n (»)B c 

= —C pT Eii (»)Ei 2 E ^2 (*)Bc + D” T C r ErJ(s)B c 
= Ci’ T E^(s)B s C c E^(s)B c + D” r C c E^(«)B e 


(23) 


The transfer function from the input noise/disturbances, <1, to performance 
output, y pr , is given by: 


W*) = GuW = r D” r C e ] 

= C’ ,r X n (»)B s + D’ ,r 
= C pT Ey l ] {.s)B, + D pr 


'Xn(s) 

A 12 (*)' 

~B a ' 

La-2! (*) 

a 22 (*)J 

_0 _ 


+ D pr 


(24) 
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• The transfer function from the exogenous disturbances, tv, to performance 
output, y pr ', is similarly given by: 

V„,( S ) = GuM = + D (25) 

• The transfer function from the tracking command, r, to control input, </, is 
given by: 

T, r (») = K(«) = C c E.J.}(s)B c (26) 

• Because the configuration is open-loop the transfer functions from the input 
noise/disturbances, <7, and , the exogenous disturbances, tv, to control input, u, 
are given by T u( { = 0 and T uw = 0. 

Now, using one of the equations (20)-(26) to calculate one of these transfer 
functions, the frequency response function matrix of the open-loop system is 
evaluated for various values of s = with u ; taking on the user-specified 
frequency values. The open-loop gain and phase plots (Bode plots) may then 
be computed directly from the frequency response function matrix, if desired. 


2.6. Efficient Calculation of Open-Loop Transfer Functions 


If these formulae are used to calculate transfer functions for a large order 
flexible system, with values of p in the hundreds or thousands, then the presence 
of the expression jET^s) effectively precludes the use of full matrix techniques, 
both for reasons of computation time and accuracy. Fortunately, the matrix 
E\ i (.s) = si - A s is block diagonal with 2x2 block diagonal elements. It is 
a common engineering practice to exploit this sparsity structure for computational 
efficiency. 

It is noted that the expression E^ (s) always occurs in one of the combina- 
tions (»)!?* or E n [ {s)B w . Additional computational efficiency comes from 
exploiting the fact that the odd numbered rows of matrices 77, and B w are zero 
(see equation (5) and the following remarks). The following illustrates how the 
sparsity is exploited. 

Assume s is not in the spectrum of A s . From equations (3) and (4), it follows 
that (si — A s )~ ] is block diagonal with the i-th block being 


( 


si - 



1 

6 -2 T 2 + Aj 


s + 




i = 1,2, p. 


(27) 
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If the row / of B s is denoted by bj \ if Q(.s) is used to represent {si - A s ) 1 B,, 
and if Q ( s ) is partitioned as 


-Qi(*r 

Q-M 


(28) 


where each partition matrix Qi(s) is a 2 x m matrix, then Q is calculated using 
the formula 


Qi( s ) — 1 1 / (*‘ 2 + ^ | 


4 2/) 


i = 1,2 p. 


(29) 


A similar approach is followed in the computation of E { [ [ (s)B tv . 

The term E.^ (*) = {si — ) 1 appears in several of the transfer function 

formulae, always in the context C'cE^ 1 (s)B r . The k x k matrix si - A c is 
generally a full matrix. Consequently, the computation of the controller transfer 
function C'cE.J.J (s)B c is achieved through conventional approaches (see, e.g., the 
work of Laub [1, 2]). 

Most of the open-loop transfer function formulae contain one of the expressions 
C 9 Eu l (»)Bs, C s E^(s)B w , 0> r E^(*)B s , or C pr E^[s)B w . By comparing the 
FLOP (FLoating point OPeration) count for alternate ways of computing one of 
these expressions, the computational advantage of utilizing the block diagonal 
structure of £ , j~ 1 1 (*) can be quantified. 

Observe that 


G 2 i(«) = = C»Q(») 

The standard, full matrix way to calculate Q(s) would involve first performing 
an LU decomposition of si - A s , followed by a backward and then forward 
solution of the triangular systems of equations using the columns of B s as right- 
hand sides. The FLOP count for this is O (//* ) + so C? 2 i(*‘) is computed 

in 0(p 3 ) +0{p 2 m ) +0(pfm) FLOPS. Thus, in the typical case where system size 
is much larger than the number of disturbances, the calculation time per frequency 
point is a cubic function of system size. 

If this calculation must be repeated for many values of s (a typical scenario), the 
technique of [1, 2] has a better average FLOP count. An initial 0(y/ 3 ) orthogonal 
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transformation must be done once; then for each -s, Q(s) is calculated in 0(p 2 m) 
FLOPs, so G -2 i(«s - ) is computed in 0{p 2 m) +0{pfm ) FLOPs. Thus, if the number 
of frequencies for which this calculation must be repeated is large enough (see 
the discussion in the last paragraph of Appendix A), the O ( /> 3 ) start-up cost is 
distributed over the frequency points such that the calculation time per frequency 
point is a quadratic function of system size. 

When the calculation is done as in equations (28) and (29), the flop count is 
0{pw), so <S'- 2 i(s) is computed in O(pfw) FLOPs. Thus, the calculation time per 
frequency point is a linear function of system size. This represents a substantial 
savings, particularly when a large number of modes is necessary to capture the 
dynamics of the system. 


2.7. Equations for Closed-Loop Configuration 

The block diagram of the system for the closed-loop configuration is shown 
in Figure 2. 



Figure 2. Block diagram for closed-loop configuration. 


The closed-loop system is more complicated. Using equations (7) and (8) and 
defining 


.4 = 


A s 

-B c Cs 


B s C c 

A c 


(30) 
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the closed-loop dynamics of the controlled structure may be written as: 




.< + 


0 B s B„. 0 

B, 0 0 -B,- 



i ; = [c» d]f;;4 

i/ r = {C pr D pr Cc]{ ' s l + [0 D> " D p ! OlJ '' 

{■I'c J I ir 


( 31 ) 


(32) 


(33) 


e = l-Cs 0]|^} + [/ 0 0 -/IT' 


■i = l 0 C c 




(34) 


(35) 


For the closed-loop system model in Eqs. (31)-(35), there are four possible inputs 
to the system (r, < 7 , u\ r) and four possible outputs (y, y pr , e, u). Therefore, a total 
of 16 transfer functions can be defined for the system, with each transfer function 
corresponding to an input/output pair. 

Irrespective of the choice of the input/output pair, the computation of closed- 
loop transfer functions from Eqs. (31)-(35) would require the inverse of matrix 
si - A at each frequency point s = jjj. Observing the closed-loop state matrix .4, 
it is obvious the block triangular form of the open-loop configuration has been 
destroyed by the coupling generated by the feedback connection of plant and 
the control system. However, the initial sparsity of the plant state matrix .4* 
is still intact. This sparsity is exploited to develop an efficient method for the 
computation of closed-loop frequency response function matrix of the controlled 
flexible structure. If sparsity is not exploited and many structural modes are 
modeled Q> is large), it follows that a large computational effort would be required 
to calculate the closed-loop frequency response function matrix, since this would 

involve the computation of the matrix term (si - . 4 ) B, where si - A is of 
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order 2 p + k. Here, B represents one of the block columns in equation (31) and the 
computation must be repeated for each value of -s = juj for all desired frequency 
values a;. 

In the following it is assumed that s is not in the spectrum of .4 (necessary for 
the transfer function even to be defined) and it is further assumed that s is not in 
the spectrum of .4.,. This further assumption is needed to enable some algebraic 
manipulation, and should not adversely affect the applicability of the following 
results. On the one hand, since .4s is the plant state matrix for a linear model of 
a flexible structure, its eigenvalues occur either at 0 (corresponding to rigid body 
modes) or in the left half plane (corresponding to damped flexible modes). On the 
other hand, it is anticipated that these results will be used to compute closed-loop 
transfer functions for s = j~ with ^ > 0. Thus, excluding the eigenvalues of .4* 
from the domain of applicability of these results does not impact the anticipated 
usage. 

The matrix term si — .4 may be written as 


'si - A» 

-B s Cc ‘ 


En(s) 

E 12 

. B c C s 

si- A c m 


. E'2\ 

E- n {s)_ 


Introduce the notation: 


A',i(») -Yi 2 (*) 
A'2l(«) A' 22 (*) 



Elli*) E\2 
E-21 E'2'2 ( S ) 


(36) 


(37) 


The assumptions which have been made about s insure that the inverses in (37) 
exist, as does Ej] 1 . Rewrite (37) as: 


A’n A* i 2 

'E it 

En' 

A 2 i X 22 

_E 21 

E’22 


(38) 


Expanding the lower left block of (38) and solving for A r 2 1 gives A 21 = 
— .Y 22 E’ 2 i£'[~| 1 . Expanding the lower right block of (38), substituting the previous 
expression for A’ 2 j, and factoring gives A 22 A = I, where A — E 2 2 — E 2 \E^ E\ 2 - 
This demonstrates that A is invertible, and justifies the application of the block 
matrix inversion formula given in [5, page 898] to find the inverse of the block 
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matrix in Equation (37): 

A = E-22 - E. n E^E VI 
An = El 1 1 + E^En^-'E-nEl' 

-Y 12 = (39) 

-Y21 = -A - 1 £- 2i £^ | 1 
A 22 = A" 1 


2.8. Closed-Loop Transfer Functions 

Formulae for the sixteen closed-loop transfer functions are derived in this 
section. Each transfer function formula is stated either in terms of the plant and 
controller system matrices and terms derived from them in the previous section, or 
as a simple alteration to a previously derived transfer function. The next section 
will address issues of computational efficiency. 

• By combining information from equations (36) and (39), it is seen that the term 
A(.s), which is used in the computation of all of these transfer functions, can 
be calculated by the formula: 


A( s ) = £ 22 ( s ) + B c C,Elt(«)B,C c 


(40) 


The transfer function from the tracking command, r, to measurement output, 
y, is given, with the aid of Eq. (37), by: 

T,r(») = (/ + G 2 l( S )AX*))-'G 2 l( S )A-( S ) 


= {C a 0] 


1 

A'i 2 (i 

O' 

1 

0 

35 

1 

X 22 (> 

OJ 

lAj 


(41) 


= C.-Yi 2 (»)B c = -C.£f l , (*)B, 2 A- | (*)B < . 

= C s £ i - 1 1 ( s )B 5 C c A-'(s)B<- 

The transfer function from the input noise/disturbance, d, to measurement 
output, y, is given by: 

Tyd(s) = (I + G 21 (a) l G' 2 i(a ) = [C s 0] 

= C S Xn(a)B s 

= C,El?(»)B, - C > Eli(»)B s C c \~ i (s)B c C > El^{a)B t 

(42) 


Yu( 4 ) .Yi 2 (*)‘ 

"B/ 

A 21 (») a 22 ( s )_ 

_ 0 _ 
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The transfer function from the exogenous disturbances, -u\ to measurement 
output, y, is given by: 


A'u(s) A| 2 (.s)’ 

’d u : 

_A'2l(s) A r 22(s)_ 

. 0 . 


T yw {s) = (I + G-2 i(*)A») G 22 (*) - [Gs 0] 

= C 8 Xn(»)B w 

= C,£r, l (*)Br - C,E 1 - 1 1 (*)B,CrA- l (*)B f C,E, 


(43) 

The transfer function from the measurement noise, r, to measurement output, 
y, is given by: 


Tyr,(s) = -Tyr(s) 


(44) 


The transfer function from the tracking command, r, to performance output, 
y pr , is given, with the aid of Eq. (37), by: 


7>>(») = G n (»)(J + G2iMA») ‘/i'M 
= [ C P r ] T-Vuts) Aiijs) 

1 |A-2i(*) A-22 (s). 

= C' ,r .Vi2(s)S< + D^CcX-aMBc 
= + D pr C c A~ , MB c 

= C pr EJ- l H»)B.C c A- l (»)B' + DrCcA-'WBc 


(45) 


The transfer function from the input noise/disturbance, d, to performance output, 
y pr , is given, with the aid of Eq. (37), by: 


7W*) = G n (*)(/ + A'(*)G 2 , (»))-’ 


= {C>’ r D pT c r 


+ D pT 


A„( s ) .Y,2(s)lp.. 

A'si(») A 22 (v)J[0 
= C pr A’i i {s)B, + D pr C c X 2l (>)B s + D pr 
= C pr E^( s )B s - crErfWB^C'A-'MBcC'ErfWB,} 

- D pr {c c A- 1 {»)B c CsE-[ l l (»)B s } + D pr 


(46) 
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• The transfer function from the exogenous disturbance, w, to performance 
output, y pr , is given, with the aid of Eq. (37), by: 

Typr w {x) = G\2(s)(I + K(s)G 2 l(s)r ' 

\r<pr DP r r' 1 12 ( • i ’) 1 

_l Crl l.Y 2 ,(») -Y 22 (s)J [ n \ +u - 

= C pT (x)B w - C'”'£'f 1 1 ( s )i? i ,{c' f A- 1 ( s )2? < C t £' 1 - 1 1 ( s )B lr } 

- D'"'{c c A- | (*)i? c C,£7 l 1 (*)S«’} + D<" 

(47) 

• The transfer function from the measurement noise, r, to performance output, 
y pr , is given by: 


Typr v (n) — —Typr r (s) 


(48) 


• The transfer function from the reference command, r, to reference error is 
given by: 

T er (s) = I - T yr (s) (49) 

• The transfer function from the input noise, d, to reference error, c, is given by: 

TM = -Ty d (s) (50) 

• The transfer function from the exogenous disturbances, w, to reference error, 
e, is given by: 


T cw (s) = - T yw {s ) 


(51) 


• The transfer function from the measurement noise, v, to reference error, e, is 
given by: 


T ev (s) = T yr ( ,) - I 


(52) 


• The transfer function from the tracking command, r, to control input, u, is 
given, with the aid of Eq. (37), by: 


T U r($) — A + G ? 2i(-‘5')A («)) — [0 C c ] 

= C c .Y 22 (s)B ( = C f A -1 (*)B c 


r-Y„(.) 

•Yi 2 (s)1 

1 

o 

£ 

i 

-Y 22 (s)J 

IAJ 


(53) 
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• The transfer function from the input noise/disturbance, d, to control input, u, 
is given, with the aid of Eq. (37), by: 

T u( i ( - s ) = — A' ( s ) ( I + G-2 i ( * ) A ( # ) ) 1 G-2 i ( « ) = [0 C ( ] 

= -C c A' 2 i(»)£* = -CcA- l (s)BcC s Eu l (s)B s 

(54) 

• The transfer function from the exogenous disturbances, w, to control input, u, 
is given, with the aid of Eq. (37), by: 


A'n(») A 12 (-s) 
A*21 (••>•) X 22 (*)J 


D s 

0 


T u w (a) = - K ( h ) ( I + G-2 i ( * ) A ( -s- ) ) 1 Gn ( - s ' ) 


! 

ri 2 (*r 

r*-i 

|x 2 i(*) : 

V 22 (-v). 

\ 

O 

i 


= -C c X 2 iG)D w = -G c A- 1 (s)D c C s E; ] 1 G)B w 


(55) 


• The transfer function from the measurement noise, v , to control input, ?/, is 
given by: 


Tu v ( * ) — T ur (s) 


(56) 


Now, using one of the equations (41)-(56) to calculate one of these transfer 
functions, the frequency response function matrix of the closed-loop system is 
evaluated for various values of $ = juj, with u; taking on the user-specified 
frequency values. The closed-loop gain and phase plots (Bode plots) may then 
be computed directly from the frequency response function matrix, if desired. 


2.9. Efficient Calculation of Closed-Loop Transfer Functions 

With the transfer functions written in this form, the following computational 
efficiencies are observed: 

□ Since Ei i = si - A s , and and B w have zeros in the odd numbered rows, the 
terms E±}B S and E^ D tr can be computed using the techniques presented for 
efficient computation of the open-loop transfer function (see the §2.6 through 
equation (29)). Achieving this efficiency was the entire point of the algebraic 
manipulations done in achieving the forms of the transfer function formulae 
given in the previous sections. 

□ The computation of A ~ X D C is done as a full matrix computation (a possible 
improvement to this is presented in Appendix A), but since A is of the same 


26 


order as the control system, which is usually small compared to the order of the 
analysis model of the plant, it should not be very costly to compute. Following 
accepted contemporary numerical analytic practice, A~ [ B ( should be computed 
by solving systems of linear equations using A as the coefficient matrix and 
the columns of B c as right hand sides ([6, Chapter 3]). 

□ The expected shapes of the matrices and the exploitation of common sub- 
expressions make it advisable not to precompute some of the matrix products 
in equations (40)-(56) which are independent of the frequency parameter s (if 
A' is a tall, skinny matrix, and Y is a short, wide matrix, so that IT' = XV 
is both tall and wide; and if c is a vector, then calculating X(Yz ) is cheaper 
than calculating TTA). 

□ Common sub-expressions, such as those mentioned in the previous items and 
those enclosed in braces in some of the transfer function equations are computed 
once per frequency, saved, and reused. Some of the intermediate results formed 
in computing A can be reused in many of the transfer function formulae. 

□ If there are no acceleration sensors in the performance outputs, so that the 
feedforward matrices D pr and D 1 ’,! are null, the software implementations of 
these computations may bypass any terms in the transfer function formulae 
which contain either of these matrices. 


The closed-loop system matrix .4 (equation (30)) has order 2p + k. As in the 
discussion of the open-loop calculation, if T y pr w (s) is calculated as 


!>„.(*) = [C« r D” r C f r ](*/-. 4 ) 


-1 


'B w 

0 


+ D pr 

' ^ W 


using standard full matrix techniques, the computation takes 0^(2y>+A*) 3 j + 

o((2 F + *) 2 «) + 0((2p + A )Uj) FLOPs per frequency point. Using the technique 
of [1] and [2], if the number of frequency points for which the calculation is to 
be repeated is on the order of 2 p + A or more, then T y pr w (a) can be calculated 

in o(^(‘2p + h fg'j + 0((2p + A )Uj) FLOPs per frequency point. Again, these are 
cubic and quadratic functions, respectively, of the system size. By counting FLOPs 
resulting from subroutine calls and DO loops in the FORTRAN software used to 
implement the calculation of equation (47), it is determined that T y pr w (s) can be 
calculated in 


0(p(lg + gf + fm + grn) + k 3 + k 2 (g + m) + k(g f + jm + gm) -f Igm) (57) 
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FLOPs per frequency point. Again, this is a linear function of the number of 
retained modes p. 

In order to evaluate the significance of Eq. (57), it is useful to recall the 
meanings of the variables in it, and to take note of their expected relative sizes 
in the applications to which this theory is intended to be applied. The number of 
retained modes p could be the largest by far of these variables, with values ranging 
even into the thousands. The number of performance outputs / and the number 
of controller states k are of moderate size, perhaps as large as multiples of ten. 
The number of disturbance inputs g is much smaller, perhaps as large as 10 or so, 
while the number of measurement outputs / and the number of control inputs to 
the plant rn are even smaller, typically being in the low single digits. Therefore, 
the terms involving j> should be looked on as the dominant part of Eq. (57). Of 
the remaining terms, the 0(A 3 ) term is the next most dominant. 

If a FLOP count expression were to be developed for another of the transfer 
function formulae, it would differ in some details. However, two important points 
would remain the same: p would appear linearly, and, of the terms which do not 
depend on />, the 0(A* 3 ) term would continue to appear and dominate. 

It is therefore interesting to note that, if the number of frequency points at 
which the transfer function is to be evaluated is large enough, the term 0(A 3 ) 
in this last expression, which comes from performing a LU decomposition on the 
matrix A, can be reduced enough to become part of the terms which are of lower 
order in k. This reduction is based on applying the technique of [1, 2] to Ep2 and 
making use of the observation that the other term in the definition of A is, in the 
expected applications, of low rank. This is seen by observing from Eq. (40) that 
the rank of the term added to £22 to get A is at most min (/, in). Details of this 
alternative calculation are given in Appendix A. 


2.10. Discrete Systems 

The computational procedure outlined for the calculation of open-loop or 
closed-loop transfer functions of the system extends to the case wherein the 
plant and the controller are represented by linear discrete-time systems with little 
adjustment. These adjustments and considerations are: 

1. The plant state and influence matrices, A s , , B w . C s , C ,pr , D pr , and Dy'\ 

which have been used in the computation of open-loop or closed-loop trans- 
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fer functions, take their discrete-time forms (instead of their continuous-time 
forms). 

2. The matrix E\ \ = si — A s continues to keep its block diagonal form with 
2x2 block diagonals. However, the special form of Eq. (27) for the inverse of 
the 2x2 block diagonal elements does not hold anymore. Instead, the general 
inverse for a 2x2 matrix is used. 

3. The influence matrices B$ and B w do not possess the property of having zeros 
in the odd numbered rows; they are full matrices with no inherent sparsity. This 
means that, in the computation corresponding to Eq. (29), both an odd numbered 
and an even numbered row of the “B” matrix are involved in the computation, 
as opposed to only an even numbered row for the continuous-time case. 

4. The frequency points s which are used in the computations are computed from 


where u/ are the desired frequencies in rad/s and T s denotes the sampling 
frequency of the discrete-time system. 

2.11. Software Implementation 

The algorithms developed in this paper has evolved over time. In the initial 
work as reported in [3], only the formulae for the open-loop and closed-loop transfer 
functions for performance output as a function of disturbance input, called T y pr w 
in this paper, were developed; and these only for the continuous case. Software 
capabilities also evolved during this time. The initial implementation of software 
to evaluate T y pr w is distinctly different from the current implementation. The 
techniques suggested in Appendix A have not yet been tried out in software. 

2.11.1. Initial Implementation 

The evaluations of the open- and closed-loop transfer function T y pr w were 
originally implemented using MATLAB function M-files (MATLAB, a product of 
The Math Works, Inc., is “a technical computing environment for high-performance 
numeric computation and visualization” [7, page i]), and as FORTRAN 77 code 
which is then accessed through MATLAB using the MEX-file external interface 
facility. The M-files have the advantage of being relatively easy to write and of 
being portable to any computer capable of running MATLAB. The FORTRAN 
MEX-files are more labor intensive to program, but for those computer systems 
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which are capable of compiling them, they perform the calculations faster. This 
is illustrated in §3. 

The M-files contain straightforward implementations of the calculations pre- 
sented in the open- and closed-loop formulae for T y vr w . The FORTRAN source 
code for the MEX-files uses the Basic Linear Algebra Subprograms (BLAS, [8-141) 
to perform vector-vector, vector-matrix, and matrix-matrix operations. In addition, 
LAPACK ([15]) subroutine zgesv, a complex double precision linear equation 
solver, is used to calculate A -1 (D c C il E'[ l l {s)B v ^ . It is in this form that the 
transfer function evaluations are implemented in version 1 of the software package 
PLATSIM [4], which is available from COSMIC (LAR- 15287). 

2.11.2. Current Implementation 

Implementing the full scope of the computational formulae developed in this 
paper was a task of significantly greater magnitude than was involved in the initial 
implementation. In the continuous case, there were seven nontrivial open-loop 
transfer functions and 16 closed-loop transfer functions for a total of 23. Extending 
the work to the discrete domain added another 23 transfer functions. Fortunately, 
there were some similarities in the forms of the calculations between some of the 
transfer functions. 

Exploiting these similarities in form, it is possible to calculate all seven 
continuous open-loop transfer functions using just four function M-files, one 
executive routine and three which do the actual calculations. Similarly, the 16 
continuous closed-loop transfer functions can be computed using just six M-files, 
one executive and five performing the actual calculations. The discrete case uses 
the same software structure, so any of the 46 distinct transfer functions can be 
calculated using a suite of 20 function M-files. 

It would have been a programming task of nontrivial magnitude to render these 
46 transfer function formulae in FORTRAN or C code from which MEX-files 
could be compiled. Fortunately, by the time this theory had been developed, The 
Math Works, Inc. had developed the MATLAB compiler [16]. This is a utility which 
translates function M-files into C code which can either be used in stand-alone 
applications (with the addition of the MATLAB C math library) or be compiled 
into MEX-files. 

There are some limitations as to which MATLAB commands can be compiled. 
Also, some calculations can be programmed in MATLAB in more than one way, 
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such that one program is more efficient as an M-file and another is more efficient 
as compiled code. So, slightly different versions of the 20 M-files are maintained 
for compilation into MEX-files as opposed to direct execution by MATLAB. These 
20 M-files are used to compile four MEX-files; one to perform the computation of 
all of the open-loop transfer functions in the continuous time case, one for closed- 
loop continuous time transfer functions, and one each for open-loop and closed- 
loop discrete time transfer functions. It is in this form that the transfer function 
evaluations are implemented in version 2 of the software package PLATSIM [17]. 
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3. Numerical Examples 


A number of numerical examples are presented to demonstrate the efficiency 
and accuracy of the algorithm presented in this paper. Execution time comparisons 
are made to two standard full matrix methods of calculating the frequency response 
function of a closed-loop system and a block elimination method for linear sys- 
tems with a multi-row and column boarder. Accuracy assessments are made by 
comparing these answers amongst themselves and with answers calculated using 
quadruple precision (128 bit) floating point computer arithmetic. Except where 
specifically noted otherwise, all software used in this study utilized ANSII stan- 
dard double precision (64 bit) floating point computer arithmetic. All tests were 
performed on a Sun SPARC 10 workstation. 

The data come from a model for the EOS-AM-1 spacecraft used in a jitter 
reduction study ([18]). This spacecraft (see Fig. 3) has a diameter of 3.5 m, a 
length of 6.8 m, and a weight of 5190 kg. 



Figure 3. Schematic diagram of EOS-AM-1 spacecraft. 
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The structural model contains 703 modes for a potential 1406 plant states. This 
modal model includes all modes with frequencies less than or equal to 200 Hz. 
The reason for using such a large flexible model of the spacecraft is to capture the 
diverse dynamic characteristics of the disturbances that act on it. The frequency of 
disturbances can vary from orbital frequency, such as gravity gradient, to sub-Hertz 
frequency, such as mass imbalance on scanning mirrors, to frequencies of more 
than 100 Hz, such as are generated by the cryocoolers. 

There are 6 rigid body modes and 697 flexible modes ranging from 1.24 to 
1564 radians per second. The 6 measurement outputs are the spacecraft’s roll, 
pitch, and yaw attitudes, and their rates, at the spacecraft navigational unit. The 
three actuators are x-, y-, and z-axis torquers. The control system has 39 states. 
Depending on which of the four types of input and four types of output were being 
used, tests were run with up to 10 channels of input and up to 27 channels of 
output. There is enough similarity among the formulae for the open- and closed- 
loop cases and between the continuous and discrete time cases that it was judged 
sufficient to conduct this testing for the closed-loop continuous time case. 


3.1. Test 1: Initial Implementation 

In this test, the closed-loop continuous time transfer function Ty P r w is calculated 
using the initial implementation of the software. Several parameters are varied 
in order to assess their effect on execution time of the software for the various 
algorithms. 

Cases were run using position measurements at each output, resulting in no 
feedforward term, and using acceleration measurements at each output, resulting 
in a feedforward term being present. 

All algorithms used in this timing study are intended to be used to calculate 
the frequency response matrix at multiple frequency points, so that the frequency 
response may be plotted (as, e.g., Bode plots). They all have some calculations 
which are done once per entry into the algorithm and other calculations which are 
done once for every frequency point. To take this into account, all cases were run 
over a range of 200 frequency points and most of them were rerun using 2000 
points (the exceptions were the cases which would have required more than 5 days 
of cpu time to complete). In all cases, the points were logarithmically distributed 
between frequencies of .01 and 10000 radians per second. 
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The executions times of test codes are compared on 12 representative problems. 
Three different plant sizes were used: a small plant with 1 input, 1 output, and 
61 states (24 structural and 39 controller states); a medium plant with 3 inputs, 
5 outputs, and 221 states (184 structural and 39 controller states); and a large 
plant with 10 inputs, 27 outputs, and 1445 states (1406 structural and 39 controller 
states). For each plant size, two different sets of outputs were used: one with no 
feedforward term (i.e., no acceleration outputs were used as performance outputs) 
and one with feedforward. The frequency response function of each of these 6 
plants was calculated at a short (200 values) vector of frequency values and at a 
long (2000 values) vector of frequency values. 

3.1.1. Software Used in Test 1 

In Test 1, two software realizations of the closed-loop frequency response 
function algorithm discussed earlier are compared to two software realizations 
of previously available algorithms for the calculation of transfer functions. 

The algorithm presented in this paper is programmed both as a MATLAB func- 
tion M-file and as FORTRAN 77 code which is then accessed through MATLAB 
using the MEX-file external interface facility. These will be called, respectively, 
the new M-code and the new MEX-code. 

One of the programs used for comparison makes use of the algorithm in [1] 
and [2]. The FORTRAN code in [2] is in single precision; Laub’s own double 
precision FORTRAN code is imbedded in the software package FREQ ([19]) and 
was used here. This test code is purely in FORTRAN 77. It will be called the 
old FORTRAN code. 

Preliminary testing indicated that, in order to get a reasonably well-conditioned 
matrix for the si — A expression^ the Laub code, it was necessary to exercise the 
built-in option of balancing the .4 matrix. The unbalanced matrix was particularly 
ill-conditioned at low frequencies. This can be attributed to the presence of the 
rigid body (0 frequency) modes. In die Laub code, balancing was coupled with 
the extraction of the eigenvalues of .4. As Laub wrote this code, the same value 
of the input flag which signaled the code to balance the .4 matrix also signaled the 
code to extract its eigenvalues. For purposes of timing tests here, the Laub code 
was modified so that the portion which extracts eigenvalues was bypassed. 

A second program used for comparison is the Math Works M-file freqrc.m, an 
undocumented utility routine in the Robust Control Toolbox, [20], which calculates 
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(to quote the on-line program preamble comments) “Continuous complex frequency 
response (MIMO)”. This will be called the old M-code. Once again, to achieve 
reasonable accuracy, it was necessary to balance A. This was done using MATLAB 
built-in routine balance. 

3.1.2. Timing Comparisons 

Particularly on the larger plants, the new algorithm performs dramatically faster 
than the older programs. This should not be taken as an indictment of the older 
algorithms. The older algorithms were designed for arbitrary plants while the new 
takes full advantage of the particular pattern of sparsity which results from using 
the modal model of a flexible structure. On the other hand, when the new algorithm 
is applicable, it enables analysis of structures of much larger order than would be 
practical or even possible before. 

Table 1 gives the time in seconds to calculate the frequency response function 
using each of the 4 test routines for each of these 12 cases (except that the old 
M-code does not attempt the two largest cases). 

One conclusion to be drawn from this table is that the timing values returned 
by the system timing software are not totally consistent with each other. The 
computations by first three software packages include an up front check to see if 
feedforward is present. The bulk of the code is executed whether feedforward is 
present or not. If feedforward is present, additional code is executed which should 
take additional time. But in 9 of 18 cases, the table shows the feedforward case 
taking less time than the one without. 

That said, there are still significant trends to be observed in this timing data. The 
new MEX-code is significantly faster than the new M-code; in the more important 
larger cases, about 3 times as fast. This justifies the effort of rendering the algorithm 
in FORTRAN and writing the interface necessary to access it through MATLAB. 

Comparing the new MEX-code, which is FORTRAN based, with the old 
FORTRAN code shows that for the small system, the old code more than holds its 
own. This is not unexpected, since in the small system, the controller dominates 
the count of states. Thus, there is relatively little sparsity from the structural part 
of the plant of which the new MEX-code may take advantage. But even in the 
medium size case, the new code is 4 to 7 times faster than the old. This is getting 
near the size at which conventional numerical analytic wisdom would place the 
limits of applicability of the old, full matrix based, technique. For large cases, 
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the new MEX-code is about 40 times faster than the old FORTRAN code. Since 
the CPU time for the old FORTRAN code is expected to grow quadratically with 
the number of system states, while that of the new technique is expected to grow 
only linearly, it is not surprising that the difference between them in the largest 
example is so great. 

In what is called here the old M-code, Math Works actually used M-files for 
outer loop logic control to drive a built in function, itifr. This function calculates 
the matrix G whose columns are (s(i)I - A)~ l h where a is a vector of complex 
numbers (set in the present application to ju>, where j 2 = - 1 and u> is a vector 
of frequencies), .4 is a system matrix (set in this application to .4), and b is a 
column vector (set in this application to a column of the D tr matrix). The on- 
line help for itifr states that it “implements, in high speed” what the user could 
calculate by looping through the elements of h and building G one column at a 
time. Despite this, it is only competitive on the smallest system, and then only 
against the new M-code which utilizes MATLAB built-in functions only at the 
more primitive level of basic matrix operations. 

All of the algorithms tested have some “once per entry” calculations in addition 
to the calculations which occur once per frequency value. Thus, the time for the 
2000 point calculations should never be more than 10 times that for the 200 point 
calculations. In Table 1, there are several exceptions to this. This reinforces the 
previous remark that the numbers returned by the computer system timing routines 
are, at best, approximate. However, from looking at the largest case, it can be 
reasonably concluded that the “once per entry” overhead is fairly small in both 
realizations of the new algorithm while being substantial in the old FORTRAN 
code, at least for large systems. This is expected, since for a system of order n 
(all other parameters being held fixed), the “once per entry” overhead in the old 
FORTRAN code includes the initial reduction of the system matrix to Hessenberg 
form which takes 0(/v 3 ) FLOPs, while the “per frequency point” calculation takes 
G(n 2 ) FLOPs. 


3.1.3. Accuracy 

No formal error analysis has been performed on the new algorithm. There is, 
however, numerical evidence to support the thesis that the new algorithm is more 
accurate than the older techniques, particularly when applied to larger systems. 
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Outputs from the four algorithm realizations were compared. For each fre- 
quency value, individual entries in the frequency response matrices computed by 
the four codes were compared using a symmetric relative error: The discrepancy 
between complex numbers // and £ (not both 0) was measured by 


Hri'X) 


\’i-Q 

•5(M + KI>' 


e symmetric relative error uses the average magnitude of the two numbers being 
compared. It is intended for use in comparing two numbers neither of which was 
known to be correct. 

This error measure ranges from a minimum of 0 to a maximum of 2. A 
value of />(//, 0 near 10“ 77 indicates that // and ( agree to about n decimal places 
while <*>(//. () > .1 indicates anything from rough approximation (near .1) to no 
correlation (bigger than, say, 1). For each fixed frequency, the worst discrepancy 
over all possible input-output pairs was observed. 

The size of the discrepancy between the frequency response function matrices 
computed by these codes was observed to depend not only on which two of the 
codes were being compared but also on the size of the system, the frequency, and 
whether or not feedforward was present. It would take too much space to present 
details of these comparisons. However, some general statements can be made. 

For each of the test problems (corresponding to one row of Table 1), the results 
produced by the four codes were compared two by two. The overall best agreement 
between any pair of calculations came from comparing the outputs of the new 
MEX-code and the new M-code. At worst, these agree to within 7 decimal places. 
This generally improves with reduction of system size or increase in frequency so 
that best agreement is within machine accuracy. No other pairing, either between 
one of the new codes and one of the old or between the two old, ever showed 
noticeably better agreement, and in general the agreement was much worse. Often, 
results from comparing the two new codes showed that the agreement of their 
computations was better than that of any other pairing by at least 2 decimal places. 
In the largest system, this advantage could increase to 5 decimal places, particularly 
for small frequencies or when no feedforward was present. 

Thus, two dissimilar implementations of the new algorithm produce results in 
good agreement. When a parallel process is applied to the older algorithm, the two 
dissimilar implementations produce results which are not in such good agreement, 
either with each other or with those of the new algorithm. 
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To provide further evidence that the results of the new algorithm are more 
accurate, the old FORTRAN code was translated to quadruple precision from its 
native double precision and was run (at a time penalty of about 32 x) on the medium 
sized problem using no feedforward and 200 frequency points. The output from 
this showed the same degree of agreement with the output from the new codes as 
they showed with each other. 

These results combine to indicate strongly that the new algorithm provides 
more accurate results than those previously available. There are theoretical grounds 
for expecting this. The^)ld way requires the solution of linear systems with the 
coefficient matrix si - A which is usually of large order. It also had conditioning 
problems which balancing ameliorated, but did not totally eliminate. 

In the new algorithm, the coefficient matrices involved in the solution of linear 
systems are A, which has the same order as the controller, and si - .4* , which is 
of order 2, for /' = 1 , • • • , p. Particularly when dealing with a large order structural 
model, the coefficient matrices used by the new algorithm are much smaller than 
the matrix si — .4 used by the old, so there is much less opportunity for round-off 
error. Any conditioning problems coming from the interaction of the frequency 
represented by s = and the i-th structural mode in the old method is isolated 
in the new method to calculating the denominator term Jf — J 2 + 2 in 
equation (29); and this only gives numerical problems if lu is so close to ujj that 
truncation occurs in forming the difference, and (i is so small that the (small) real 
part Af — to 2 is a significant part of the whole term. 

3.2. Test 2: Current Implementation 

In Test 2, seven different closed-loop continuous time transfer functions are 
calculated using the current implementation of the software. It was deemed 
unnecessary to consider as many combinations of the parameters as was done 
in Test 1. 

As noted in §2.11.2, the 16 closed-loop continuous transfer functions are 
calculated using one of five MATLAB M-files. Two of the five accept feedthrough 
matrices as an optional parameter and conditionally perform calculations which 
depend on the feedthrough matrices. The seven test cases are chosen so that all 
five M-files are exercised, and the two which accept feedthrough matrices are 
exercised both with and without the optional feedthrough matrices. 
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It was decided that the phenomenon of algorithm start-up time was adequately 
examined in Test 1. So, in Test 2, a 301 point frequency vector with the points 
distributed logarithmically between frequencies of .01 and 10000 radians per second 
was used for all cases. 

Linear systems with two different numbers of state variables were used in 
calculating each of the seven transfer functions; one had the same number of states 
as the medium plant of Test 1 with 221 states (184 structural and 39 controller) 
and the other had the same number of states as the large plant of Test 1 with 
1445 states (1406 structural and 39 controller). The numbers of inputs and outputs 
and the presence or absence of feedforward depended on the individual transfer 
function and are given in the following chart: 


Transfer function 

Number of inputs 

Number of outputs 

Feedforward 

present 

Tyr 

6 

6 

no 

T y „. 

10 

6 

no 

TyV T r 

6 

27 

yes 

TyPr w 

10 

27 

yes 

Ter 

6 

6 

no 

Tur 

6 

3 

no 

T U cI 

3 

3 

no 


3.2.1. Software Used in Test 2 

In this study, two software realizations of the closed-loop continuous time 
frequency response function calculation algorithm given in this paper are compared 
to two software realizations of previously available frequency response function 
calculation algorithms and to one special purpose linear equation solver. 

As noted in §2.11.2, the current software implementation of the calculation of 
the 16 closed-loop continuous time transfer functions is programmed as a collection 
of six MATLAB function M-files. These will be referred to as the new M-code. 
Minor modifications of these M-files are then translated into computer code in 
the C programming language using the MATLAB Compiler. The C code is then 
compiled into a MATLAB MEX-file using the MATLAB script cmex. This will 
be called the new MEX-code. 
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The same old FORTRAN code and old M-code used in Test 1 are used here 
for timing comparisons. Also, a special purpose linear equation solver due to 
W. Govaerts and J. D. Pryce [21] is used for both timing comparisons and accuracy 
assessment. 


The Govaerts-Pryce (GP) algorithm assumes that the user has a matrix and 
further that the user knows how to solve linear systems using both that matrix 
and its transpose as the coefficient matrix. GP further assumes that what the user 
really wants is to solve linear systems using a larger matrix formed by adding some 
rows and columns to the given matrix. By using information from user provided 
solutions to linear systems with the smaller matrix and its transpose as coefficient 
matrices, GP calculates the solution to the larger system. 


In Test 2, GP is used to calculate expressions of the form ( si — 



-1 _ 

D (see 


§2.7). The block structure of the coefficient matrix si— A shown in equation (36) is 
used. It is easy to solve equations using the upper left comer E[ j (Y) or its transpose 
as a coefficient matrix since E n(Y) is block diagonal with 2x2 blocks. Since 
the GP calculation is being included in the comparisons primarily as a baseline 
for accuracy checking, equations with these 2x2 blocks or their transposes as 
coefficient matrices are solved using full pivoting Gaussian elimination, which is 
at least as accurate as the Cramer’s rule based solution process given in equation 
(27) for use in implementing the present algorithm. 


3.2.2. Timing Comparisons 

First, take note of comparisons which will not be made. Tests 1 and 2 were 
made about 2 years apart. While they were made on the same computer, there were 
subsequent hardware and software upgrades to the computer. This means that, 
even if the tests had been identical, the results would not have been comparable. 
Furthermore, while the same data set is used to generate test cases for both tests, 
none of the 12 cases used in Test 1 exactly matches any of the 14 cases used in 
Test 2. So, no comparison will be made between actual running times of any of 
the Test 1 cases with those of any of the Test 2 cases. 

Evaluating seven different transfer functions using systems of two different 
sizes for each gives a total of 14 test cases. These cases were evaluated using the 
software described in the previous section with the exception that the old M-code 
was not applied to the large systems since this would have taken an unreasonable 
amount of computer time. The cases using M-code or MEX-code were timed 
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using the MATLAB built in function cputime. The old FORTRAN code and the 
GP code were timed using Sun FORTRAN timing function dtime. Results of the 
63 individual timing runs are given in Table 2. 

The current implementation of the algorithm continues to outperform traditional 
software, with the degree of improvement increasing with the size of the plant. On 
the medium sized test cases, the new M-code ran faster than the old M-code by a 
median factor of more than 34, while the median speed up of the new MEX-code 
over the old FORTRAN code was by a factor of more than 2.3. No attempt was 
made to run the old M-code on the large cases. On the seven large cases, the new 
MEX-code showed speed ups over the old FORTRAN code by factors ranging 
from 42 to 235 (median, 153). 

The new MEX-code still outperforms the new M-code. In the medium sized 
cases, the median speed up of new MEX-code over new M-code is by a factor of 
1.5 and in the large sized cases, the median speed up factor is 1.75. This speed 
up is not as good as in Test 1. 

The difference in speed up factors is probably attributable to the way in which 
the respective MEX-codes were generated. In Test 1, the MEX-code was compiled 
from handwritten FORTRAN code. In Test 2, the MEX-code was compiled from C 
code which was generated from M-files using the comparatively recently released 
MATLAB utility mcc. So, the Test 1 code had a couple of advantages. It was 
hand written, which has at least the potential of producing more efficient code than 
an immature machine translation scheme. It was written in FORTRAN, and, in 
the authors’ experience, contemporary FORTRAN compilers seem to do a better 
job of optimizing scientific calculation code than do contemporary C compilers. 
The decision to use the automatically generated code instead of hand-crafted code 
for the MEX-files in the current software implementation of the algorithm of this 
paper was an economic one based on available man-hours. 

Although testing with the Govaerts-Pryce code was included primarily to 
assess the accuracy of the current implementation of the new algorithm, timing 
measurements were taken. The new MEX-code was also faster than the GP code, 
with median speed up factors of 8.9 for the medium sized plants and 40.6 for the 
large plants. So, while the GP code takes advantage of the sparsity in the structural 
plant matrix in its own way, the algorithm of this paper is still preferable for 
computational speed. While the old FORTRAN code was designed for efficiency 
in computing (full matrix) transfer functions at multiple frequency points, on the 
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large plant, the GP code, with its exploitation of sparsity, was faster than the old 
FORTRAN code by a median speed up factor of 3.78. 

3.2.3. Accuracy 

Numerical evidence from the results of Test 2 continue to support the thesis that 
the new algorithm is more accurate, particularly on large problems, than traditional 
transfer function calculation techniques. When Test 2 was made, the availability of 
the GP software not only gave another independent calculation of transfer function 
matrices for use in comparison, but also gave a new tool for generating a measure 
of relative error. 

As explained in Appendix B, the discrepancy between the results of calculating 
a transfer function 

T = c(sI-Xy X B + D 

in two different ways is measured on a component by component basis relative to 
the following quantity which will be called a comparison magnitude: 

1 

The challenging part of computing this comparison magnitude, the determination of 

{^1 ~ ^ > was accomplished using GP to solve the system of linear equations 

having a I — A as a coefficient matrix and the identity matrix on the right hand 
side of the equation. All values for discrepancies and errors reported subsequently 
are relative to this comparison magnitude. 

As in Test 1, unless otherwise mentioned, all computations were performed in 
64 bit ANSI standard floating point computer arithmetic, called “double precision” 
on the Sun workstation used for this study. In Sun double precision arithmetic, 
machine epsilon, the smallest number which, when added to 1 in computer 
arithmetic, produces a result which is different from 1, is 2~ 52 which is about 
2.2x 10~ 16 . 

In order to judge the accuracy of the GP code on the test problems, a quadruple 
precision (128 bit) version was created and exercised on the same test problems 
(using every 10 th frequency point; the median over the 14 test cases of computation 
time per frequency point for quadruple precision was about 70 times that for double 
precision). Except for the two test cases calculating T ur , the discrepancy between 
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relative error 


the double and quadruple precision results was never more than 10 , and was 

usually less than 10“ 15 . Even in the exceptional cases, the discrepancy never 
exceeded 10“ 12 and was usually less than 10 -14 . Because of these comparisons, 
the GP results were taken as the accuracy standard against which other results 
were measured. 
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Figure 4. Discrepancy between the transfer function values calculated 
by the new MEX-code and the GP code, all test cases shown. 

Of particular interest in this paper are results reflecting on the accuracy of the 
algorithms developed herein. The two software realizations of that algorithm, the 
new MEX-code and the new M-code, can be viewed as being compilations by two 
different compilers of essentially the same code. It is expected that they would 
give very similar results. Indeed, when compared, discrepancies between their 
results were uniformly bounded over the 14 test cases by 10 ^ with a substantial 
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proportion being identical. So, it is enough to judge the accuracy of the present 
algorithm to compare the results computed by the new MEX-code to the results 
computed by the GP code. 

The 14 test cases each calculated transfer function matrices at 301 different 
values of frequency. These matrices had dimensions ranging from 3x3 to 10x27, 
depending on the case. A total of 355,782 data points were generated. A scatter 
plot showing the discrepancies between the values calculated by the new MEX- 
code and those calculated by the GP code at all 355,782 data points is given in Fig. 
4. Of those, 88.69% of the new MEX-code answers showed an error (compared 
to the GP answers) of less than 10 -14 , 99.65% showed an error of less than 10 -11 , 
99.984% showed an error of less than 10~ 8 , and 99.987% showed an error of less 
than 10 -6 . The remaining 46 data points (.013%) showed errors of at least .01, with 
the worst case discrepancy being 2.8xl0 17 . These errors were troubling enough 
to warrant further investigation. Results of this investigation are reported at the 
end of this section. 



Figure 5. Transfer function matrix element T pr ( 13,2) for the 
large plant calculated by the new MEX-code or the GP code. 

The transfer function matrix element T pr { 13, 2) for the plant with 703 structural 
modes is a typical result of these calculations. A Bode plot of this transfer function 
as calculated by the new MEX-code or the GP code is shown in Fig. 5. The same 
plot can be shown for both calculations since the discrepancy between them, as 
shown in Fig. 6, is below the plot resolution. 
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Figure 6. Discrepancy between the transfer function matrix element 13. 2) 
for the large plant calculated by the GP code and the new MEX-code. 
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Figure 7. Transfer function matrix element T pr { 13,2) for 
the large plant calculated by the old FORTRAN code. 

With respect to accuracy, the old FORTRAN code and the old M-code showed 
very similar accuracy behavior, so only the old FORTRAN code will be discussed. 
For two of the transfer functions, T yir and T pw , accuracy of the old FORTRAN 
code was adequate, with errors of 10 -7 or better. However, for the remainder of 
the transfer functions, the accuracy results were frequency dependent. Accuracy 
was good to adequate for frequencies between .01 and 1 rad/s. It then tended to 
degenerate as frequency increased from 1 to 10, then ramped up sharply so that 
when frequency reached 10 4 rad/s, the size of errors ranged from 10 8 to 10 36 . 
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Observe the same transfer function matrix element T pr { 13, 2) for the same plant 
which was considered in connection with calculations using the GP and MEX- 
codes. A Bode plot of this transfer function as calculated by the old FORTRAN 
code is shown in Fig. 7. Over the last decade of frequency values, this plot diverges 
noticeably from the plot in Fig. 5. The magnitude of the discrepancy between these 
two calculations is shown in the scatter plot in Fig. 8. 



frequency, rad/s 


Figure 8. Discrepancy between the transfer function matrix element T pr ( 13, 2) 

for the large plant calculated by the GP code and the old FORTRAN code. 

The exceptionally large errors. 

Recall that a small fraction of the new MEX-code calculated data points, 46 
out of 355,782 or .013%, showed errors of at least 0.01 when compared to the 
GP code calculations. The worst case discrepancy was 2.8xl0 17 . This point was 
truly exceptional. The next worst discrepancy was less than 1.6xl0 2 , and only a 
total of 7 errors were more than 1. 

These larger errors occurred in several different transfer function evaluations 
(Typr r , T rr , T ur , and T U( [), but all in cases involving the 703 mode structural model, 
and all at the same frequency, 794.3282 rad/s. The worst error occurred in the (1,4) 
position of the transfer function matrix T ur , so it was chosen for further analysis. 
Recall, Eq. (53), that 

T*M = C' A-'( S )B ( .. 

The baseline Bode plot for the transfer function matrix element T„ r (l,4) for 
the large plant as calculated by the GP code is shown in Fig. 9. When the same 
calculation is done by the new MEX-code, it results in the Bode plot shown in Fig. 
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Figure 9. Transfer function matrix element T ur (l. 4) 
for the large plant calculated by the GP code. 



Figure 10. Transfer function matrix element T ur (l,4) 
for the large plant calculated by the new MEX-code. 


10. The numerical problem shows up as a spike on the magnitude plot at about 800 
rad/s. It can also be seen in the scatter plot of discrepancies between these two 
calculations shown in Fig. 11, where the exceptionally bad point is emphasized 

with an X. 

Again, this is the exceptionally bad case. The spike on the next two worst 
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frequency, rad/s 


Figure 11. Discrepancy between the transfer function matrix element TV,. (1,4) 
for the large plant calculated by the GP code and the new MEX-code. 

cases would be about 1/8 as tall as this one, with the next worse after that being 
about 1/16 as tall. 

The calculation of TV (1,4) was repeated using the new MEX-code and the 
GP code and a grid of 30,001 points; i.e., 100 times the sampling density of the 
previous calculation. The error occurred at other points of this finer grid, but 
all close to the point already discovered. A detail plot of this area covering all 
additional erroneous points is shown in Fig. 12. The points from the magnitude 
plot of Fig. 10 are marked here with an x. The baseline GP calculation is plotted 
using a solid line, the new MEX-code calculation using a dotted line. Additional 
error points were found; the magnitudes of errors were all about the same and the 
frequencies at which they occur were limited to the range from about 750 rad/s 
to 810 rad/s. 

The actual calculation of T ur (s) was being done by a MATLAB statement of 
the form: 

>> T = Cc* (DELTA\Bc) ; 

Here, the MATLAB variables cc and Be hold the values of the matrices C c and 
B e , respectively, while the MATLAB variable delta holds the value of the matrix 
A(.s) with s = 794.3282V The (1,4) element of T ur can be considered as being 
the result of the MATLAB calculation 

>> T (1, 4) = c*(DELTA\b); 

where c is row 1 of cc and b is column 4 of bc. A wide range of magnitudes was 
observed in this calculation. The calculated values in the column vector deltas 
varied over about 47 orders of magnitude. When the product c* (DELTA\b) was 
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Figure 12. Detail of fine grid calculation of transfer 
function matrix element T ur ( 1,4) for the large plant. 


formed, zeros in the row vector c canceled out elements of the top 16 orders 
of magnitude in deltas. Thus, for t(4,d to be accurate, small components of 
DELTA\b needed to be computed accurately, and in the case under study, this did 
not happen. 

Several experiments were conducted to see if improving the linear equation 
solution process implicit in the MATLAB computation delta\bc would give values 
of T ur which agree with the GP solution. Two were successful at all but one data 
point; that one data point being the same (1,4) element of T ur at frequency 794.3282 
rad/s considered previously. When two cycles of iterative improvement were used 
to improve the accuracy of the computed value for A - 1 B c , except for the one 
point noted, the results came into agreement with the results produced by the GP 
code. The same agreement was achieved when A ~ l D c was computed using the 
linear equation solution subroutine zgesvx from LAPACK [15], Even for the one 
bad point, the error decreased from 2.8 x 10 17 to 2.0x 10 2 . So, the problem of the 
few inaccurate points was traced to numerical behavior of the MATLAB linear 
equation solver and could generally be relieved by using a more robust linear 
equation solver to calculate expressions involving A -1 . 

When the complete suite of tests was recalculated with the new MEX-code 
modified to use LAPACK subroutine zgesvx for all computations involving A -1 , 
the discrepancies of the results, when compared to those of the GP code, are shown 
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Figure 13. Discrepancy between the transfer function values calculated by the 
new MEX-code with ZGESVX and the GP code, all test cases shown. 

in Fig. 13. By comparing this figure to Fig. 4, it can be seen that not only has 
the worst point improved and the rest of the violators of the 10 -2 relative error 
magnitude line been reduced to levels consistent with the GP results, but also there 
has been a general improvement in accuracy in the lower frequencies. The down 
side to the accuracy improvement using zgesvx is that the median time over the 
14 cases to calculate the transfer function matrices using zgesvx is 1.7 times that 
of using the MATLAB matrix division operator. 

It was also observed that this phenomenon occurred at a frequency where the 
affected transfer function values had experienced a large amount of roll-off. From 
a point of view of engineering analysis, even if these inaccuracies remained in 
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the calculation, they would probably have no impact on any engineering decisions 
made on the basis of this calculation. 
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4. Concluding Remarks 

An efficient and novel procedure has been developed for the calculation of the 
frequency response function of a large order, flexible system implemented with a 
linear, time invariant control system. The procedure takes advantage of the highly 
structured sparsity of the system matrices of the plant in normal mode coordinates. 
This reduces the computational cost from a quadratic function of the order of the 
system to a linear function, thereby permitting the practical frequency analysis of 
systems of much larger order than by traditional, full-matrix means. Formulations 
have been given for both open and closed-loop systems for both continuous time 
and discrete time formulations. 

Numerical examples were presented wherein the advantages of the present 
algorithm over traditional approaches, both in speed and in accuracy, have been 
demonstrated. When the initial software implementation of the algorithm was 
exercised on the largest systems, the MATLAB MEX-file based on a FORTRAN 
implementation of the new algorithm (the new MEX-code) was about 40 times 
faster than the old FORTRAN code when many frequency points were used while 
the advantage increased to a factor of about 80 or better when the calculation 
involved few frequency points. In this latter case, MATLAB M-files based on 
the new algorithm (the new M-code) was over 200 times faster than the old 
M-code. On the largest systems, the current MEX-code implementation of the 
algorithm, C language code generated automatically from MATLAB M-files by 
the MathWorks mcc, was faster than the old FORTRAN code by factors ranging 
from 42 to 235 depending on which input/output pair and transfer function matrix 
size were chosen. Even though the Govaerts-Pryce (GP) code, which was used 
for an accuracy baseline, took advantage of the same sparsity in the plant system 
matrix as did the new MEX-code, the new MEX-code was faster than the GP code 
by a median factor of 8.9 for medium sized plants and by a median factor of 40.6 
for large sized plants. 

In this study, the standards of accuracy were set by the GP code. With only 
a tiny number of exceptions (0.013% of all data points), results calculated by 
the new MEX-code and the new M-code agreed with the GP code results to at 
least 8 decimal places. The exceptionally large errors were traced to problems 
experienced by the MATLAB linear equation solver with certain of the A matrices 
for coefficient matrices. By contrast, at a substantial fraction of data points, the 
results of the old FORTRAN code and the old M-code disagreed with those of the 
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GP code by totally unacceptable error amounts, with discrepancies ranging up to 
many orders of magnitude. 
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Appendix A: Faster Calculation of C c A~ \s)B c 

The expression C c A~ l (s)D c occurs (directly or indirectly) in the formulae for 
all 16 of the closed loop transfer functions presented in this paper. If standard full 
matrix techniques are used to calculate this expression (the product A _ 1 ($)Z? C is 
calculated by solving a linear system with A(.v) as the coefficient matrix and B c 
as the right hand side and the calculation is completed by a matrix multiplication), 
then the FLOP count for this calculation is Q(k l + A' 2 / + Am/). This appendix 
presents a technique which can be used to replace the A - 3 term by A 2 m if transfer 
functions must be calculated at enough (order of A / ( / + m ) or more) different 
values of s. 

First, some mathematical machinery must be developed; then it is applied to 
this specific problem. 

A.l. Solving Low Rank Modifications of Linear Systems 

In this section, it is shown how to express the solution of one system of linear 
equations in terms of the solution of another system of linear equations whose 
coefficient matrix is a low rank perturbation of the coefficient matrix of the original 
system. The conventions which have been established in this paper for the use of 
symbols are suspended for the remainder of this section in favor of: 

Notation: 

I Identity matrix of size appropriate to context 

All other capital letters Matrices over some field 

x and y vectors over some field 

All other lower case letters Positive integers 

The results of this section are valid for matrices with entries from an arbitrary 
field. For engineering purposes, the most important fields are the real and complex 
number fields. 

In this section, A represents an arbitrary non-singular matrix. The results of 
this section are most useful when *4 has the property that linear systems using .4 
for a coefficient matrix can be solved more easily than would be the case for an 
arbitrary matrix of the same size as .4. For example, .4 could have a known LU 
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decomposition (factorization into lower triangular and upper triangular matrices), 
A could be block diagonal, or (the case of interest in this paper), .4 could be upper 
Hessenberg (i.e., -4 is zero below the first subdiagonal). 

The main result of this section depends on a lemma. 

Lemma: Let A be a non-singular k x k matrix, and let R and S be A- x r matrices. 
Then, -4 + RS T is non-singular if and only if I A S T A~ l R is non-singular. 

Proof: Both parts of the “if and only if’ will be proven by contradiction. 
Suppose that .4 + RS T is singular. Then, there exists an x ± 0 for which 
(A + RS t )x = 0. Then, .4./: = R(-S t x) . Since .4 is non-singular, At ± 0, and 
so -S T x ± 0. If I A S T A~ l R is applied to this nonzero vector, the result is 0: 

[i + S t A-'e)(-S t .) = -S T x + S T A-'A. r 

= —S T x A S T x = 0 

This demonstrates that / 4- S ^ .4 ' R is singular. This establishes the if portion 
of the Lemma. 

Now, suppose that I-\-S^ A~ * R is singular. Then there exists a ij 0 for which 
(I + S t A~ [ R)ij = 0. Then, y = S T (-A~ l Ry). Therefore, -A~ l Ry ± 0. If 
A A RS t is applied to this nonzero vector, the result is 0: 

(.4 + 7?S T )(-^- 1 /f«) = -Ru + ES T ( y -A~ l Ii!/) 

= —Ry + Ry = o 

This demonstrates that *4 + RS^ is singular. This establishes the only if portion 
of the Lemma. □ 

The Low Rank Modification Theorem for Solution of Linear Equations , which 
is also a corollary to the Sherman-Morrison-Woodbury formula [6, p. 51] or [22, 
124], can now be stated: 

Theorem: Let .4 be a k x k matrix, R and S be k x v matrices, and D be a 
A x m matrix. Suppose that .4 and .4 + RS f are non-singular. Then there exist a 
k x r matrix Z and k x m matrices Y and W satisfying, respectively, the equations 
AZ ~ R, AY = D, and (I A S T Z)W = S T Y. Further, if A' = Y - ZW, then 

A satisfies the equation (A A RS T )X = B. 

Remark: Note that AaRS t is perturbed from .4 by a matrix of rank at most r. 
This theorem says that one can solve m linear equations with the order k coefficient 
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matrix .4 -f RS 1 by solving in + r linear equations with the order A coefficient 
matrix .4 and m linear equations with the order r coefficient matrix / + S T Z. If 
the matrix .4 has properties which allows one to solve equations with coefficient 
matrix .4 more rapidly than those with coefficient matrix .4 + RS T , and if r is 
substantially smaller than A , and if in is not too large, then it might be significantly 
more efficient to calculate X as A' = Y - Z\V than to calculate A' by solving 
(-4 + RS r )X = D directly. 

Proof: The non-singularity of .4 insures the existence of Z and Y' with the 
stated properties. Since I + S T Z = / + S^A~ 1 T, the lemma insures that I + S T Z 
is non-singular, so there exists a TF with the stated property. Let X = Y — ZTF. 
It follows from the equation defining IF that S T ZW = S J Y - TF. Therefore, 

(-4 + RS T ^j X = (.4 + i?S r ) (Y - ZW) 

= AY- AZW + RS t Y - RS t ZW 
= D - i?TF + RS t Y - r(s t Y TF) 

= B 

□ 


A.2. Application to Calculation of C c A 1 (s)B c 

The idea presented by Alan Laub in [1, 2], as applied to calculating the transfer 
function of the control system, would replace the control system equations (8) by 

■i'h ~ Aj/Xh + Bj t e 
u = C h x h 
e =r - y - v 

where Ah = T~ 1 A C T, B/ t — T~ [ B C , and C), = C c T\ and T is an orthonormal 
matrix chosen so that the matrix Ah is upper Hessenberg (has zeros below the first 
subdiagonal). Essentially, this just performs a change of basis in the controller 
state space without changing the controller itself. The computational burden of 
performing this transformation is 0(A 3 ) FLOPs, but it needs to be done only once 
even if transfer functions are to be evaluated for multiple values of s. 

If this transformation is not done, then to calculate the transfer function of 
the controller, it is necessary to evaluate (si - A c )~ l B c for each s. Computing 
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this expression requires 0 (fr 3 + £ 2 /) FLOPs. If, however, the transformation 
is performed, then, along with . 4 /,, si - A/, is also upper Hessenberg, and the 
corresponding evaluation of (si - Aj t )~ l B c requires only 0 (k 2 f) FLOPs. 

It was already suggested in the discussion of computing open-loop transfer 
functions that this technique be applied to computing the open-loop expression 

C c X 22 ( s )B c = CcEn\s)B c = Cc(sl - Ac)~ 1 B c = Cl, (si - A h )~ l B k 

Operationally, the simplest way to accomplish that is to replace the original 
representation of the control system by its equivalent representation using an upper 
Hessenberg “A” matrix. 

The problem with calculating closed-loop transfer functions is that the open- 
loop relation A’22 = £22 * s replaced, in the closed-loop case, by the relation 
A22 = A -1 . Now, recall from Eq. ( 40 ) that 

A(») = EtiU) + B c C,Ert(s)B,C' 


Furthermore the rank of the term added to £22 to get A is at most min ( /, in). 
It is now a small step to write A = £22 + BS 1 where either B = B c and 
S T = C s E; l l (s)B s C e or B — B c C s E^(s)B s and S — C c depending on whether 
/ or w were smaller, respectively; and apply the Low Rank Modification Theorem 
for Solution of Linear Equations to the calculation of A ~ 1 (s)B c . If A c is upper 
Hessenberg, so that the solution of linear equations with coefficient matrix £22 (•'>') 
can be calculated substantially faster than for arbitrary k x k matrices, then the 
calculation of A - 1 ( s)B c will be accomplished with fewer FLOPs than if full matrix 
techniques were used. The calculation of C c A ~ 1 (a)B c is now completed with a 
matrix multiplication. 


A.3. Comparative Computational Complexity 

In this section, the computational burden of computing C c A~ l (.s)B c by stan- 
dard full matrix techniques is compared to the computational burden of computing 
the same expression using the low rank modified Hessenberg technique of the 
previous section. Specifically, FLOP count expressions will be developed for the 
amount of calculation required to pass from the point at which C 6 E U [ {s)B s has 
been calculated to the point at which the calculation of (s)B c is completed. 
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Full matrix FLOP count 

In this FLOP count, it is assumed that A c and, therefore, £22 (s) = si - A c 
are full matrices. 

The term B c C s E^(s)B s C c can be calculated from C ii E^(s)B ii in 
0{kfm + k 2 m) FLOPs by doing the D c multiplication first. The addition 
of E -22 to form A does not add significantly to the FLOP count. Calculating 
a partially pivoted LU decomposition of A requires 0(A 3 ) FLOPs. Using this 
LU decomposition, the calculation of X~ l (.s)B c can be completed in 0(k’ 2 f) 
FLOPs. The matrix multiplication required to calculate C c A~ [ (s)B c requires 
O(kfm) FLOPs. Thus, the total FLOP count to calculate C c A ~ 1 {s)B r from 
C s E^{s)B a using full matrix methods is 0(A 3 + k 2 (f + m) + A/m). In the 
expected applications of this algorithm, k is significantly larger than / or in, so 
the dominant term of this FLOP count formula is 0(A- 3 ). 

Low rank modified Hessenberg FLOP count 

In this FLOP count, it is assumed that A c and, therefore, £22 (•'>') = si - A c are 
upper Hessenberg matrices. The FLOP count will be based on applying the Low 
Rank Modification Theorem for Solution of Linear Equations (LRMTSLE). These 
expressions will be developed under the assumption that m < f. If f < in, the 
alternate expressions given earlier for R and S are used, and the roles of rn and 
/ are reversed in the FLOP count expression, but this does not alter the dominant 
terms. 

In order to express A in the form £ 2 2 + BS T needed to apply LRMTSLE, R is 
taken to be the k x rn matrix £ r C. s £ n 1 B v at a computational cost of O(kf rn), and 
S T is taken to be C c at no computational cost. Because £>2 is upper Hessenberg, 
the (now A x m) matrix Z of LRMTSLE can be calculated as the solution to 
the equation £ 22 Z = R in 0{k 2 rn ) FLOPs. The (now A x /) matrix Y of 
LRMTSLE can be calculated as the solution to the equation E 22 Y = B c in 
0(k 2 f) FLOPs. The matrices S^Z and S^Y needed to calculate TT can be 
calculated in 0{km 2 + km f) FLOPs. The calculation of the rn x / matrix TF as 
the solution of the equation (7 + S T Z)W = S T V must be done by full matrix 
techniques. This requires 0(m 3 + f m 2 ) FLOPs. The matrix X of LRMTSLE, 
which is A ~ l (s)B c , is calculated as X = Y—ZW in O(kf m) FLOPs. The matrix 
multiplication required to calculate C c X~ l {s)B c requires O(kfrn) FLOPs. Thus, 
the total FLOP count to calculate C c A~^(s)B c from C s E^(s)B b using the low 
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rank modified Hessenberg method is 0(k 2 (m + /) + k(fm + rri 2 ) + fin 2 + m*). 
In the expected applications of this technology, k is significantly larger than / or 
rn, so the dominant term of this FLOP count formula is 0(k' 2 (m + /)). With 
the expected values of the parameters, this is a significant improvement over the 
full matrix formula. 

This still leaves the question of how the total computation time is affected by 
adding in the time necessary to reduce an arbitrary controller representation to one 
in which the system matrix is in upper Hessenberg form. The general answer is 
that if there are enough distinct frequency points over which to amortize the cost of 
the Hessenberg reduction, then adding the cost of this initial reduction to the cost 
of computing the transfer function for all frequency values still results in a cheaper 
calculation than using the full matrix method. But, that brings up the question of 
how many frequency points are enough. One answer to this is given by observing 
that if there are A /( / + rn) or more points, amortizing the dominant 0(1" ) term 
of the start-up cost over these points contributes an additional 0(0 Hf + m)) or 
less to each. This absorbs into the existing dominant term. Another answer can be 
given by using explicit FLOP counts given in [1]. The dominant term in the FLOP 
count for the orthogonal Hessenberg reduction (the more numerically stable of two 
techniques mentioned in [1]) is 5A' 3 /3. The dominant term in the FLOP count for 
the full matrix solution method is A */3. This suggests that the break-even point for 
using the low rank modified Hessenberg method is that the calculation needs to be 
done for at least five different values of a (i.e. at least five different frequencies). 
In typical applications (such as making a Bode plot), there are many more than 
five frequencies at which the transfer function is evaluated, so using the low rank 
modified Hessenberg method should give the answers with fewer FLOPs. 
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Appendix B: Numerical Error Measurement: 

The Data Relative Error 

Numerical calculations are subject to many kinds of error. Typically, only a 
fixed precision, a fixed number of decimal places or binary bits, are carried through 
the calculation, so round-off becomes a source of error. If two nearly equal numbers 
of a given precision are subtracted, the cancellation of high order digits or bits 
results in an answer of lower precision. An algorithm which would, theoretically, 
require infinitely many steps to arrive at the correct answer is terminated after only 
a finite number of steps producing a truncation error. In a long calculation, such 
errors can be cumulative, and can be magnified by ill conditioning of the problem. 

In order to experimentally estimate the magnitude of such error in some given 
calculation, one might find oneself comparing two numbers to see if they are close 
together. For example, if one has performed a computation which is supposed 
to have solved an equation, one might substitute the computed candidate for a 
solution back into the equation, and see how close the left and right sides are. Or, 
one might use two (or more) different algorithms to approximate the same result 
and compare these approximations with each other. The question then is whether 
two numbers are “close” to each other. 

There is no absolute answer to this question. For example, the absolute error 
between 15,995 and 15,999, i.e., the magnitude of their difference, is more than 
the absolute error between .001 1 and .0023. But most people would judge that 
the numbers of the first pair are almost the same while the numbers of the second 
pair are substantially different. This is based on looking at the relative error, the 
absolute difference of the numbers relative to some appropriate magnitude such as 
their sizes. The symmetric relative error measure, <5, introduced in §3.1.3, reflects 
this judgement, since <5(15995, 15999) « 0.00025 while £(.0011, .0023) « 0.71. 

The problem in calculating a relative error between two numbers becomes one 
of finding the comparison magnitude, a quantity of an appropriate magnitude to 
which the absolute error may be compared. If one of the numbers is known to be 
the correct one, and it is not zero, then its magnitude may be used; this is what is 
commonly known as relative error. The symmetric relative error uses the average 
magnitude of the two numbers being compared. It is intended for use in comparing 
two numbers neither of which was known to be correct. 
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However, f> is not always a good measure of whether two numbers which are 
supposed to be approximately the same really are. For example, suppose .< « is an 
approximate (e.g., calculated) solution to the simple scalar linear equation ax = b 
which can also be written ax-b = 0. If .r„ is substituted into either of the equations 
for which it is supposed to be a solution, one expects the left and right sides to 
be close to each other. Suppose that the result of the machine calculation ax u is 
not exactly b, and suppose that b is used to measure this closeness. Using the first 
equation, the closeness will be measured as b{ax u J>), which will probably be on the 
order of machine epsilon. However, using the second equation, the closeness will 
be measured as 6(ax a - M) = 2, which is the maximum measure of discrepancy 
which the symmetric relative error measure can return. This second test has given 
a false sense of the quality of x a as a solution to the equation ax - b = 0. 

A problem with trying to measure how good an approximation x a is as a 
solution to the equation ax — b = 0 using b(ax n — b, 0) is that this symmetric 
relative error calculation does not take into account the magnitude of the data 
which went into calculating x„. One way to take data size into account would be 
to measure the error relative to the size of the constant term b. As long as b / 0, the 
relative error measure becomes \ax u — b\/\b[, the error is being measured relative 
to the magnitude of the constant term of the equation. 

This actually gives a good indication of how close the residual ax n — b is to 
0. However, if the scalar equation ax — b is replaced by a system of equations 
represented in matrix-vector form as Ax = b, the natural extension of this relative 
error measure may not give such a good idea of whether all components are near 
0. This extension would be to take \\Ax a - b\\/\\b\\ as a measure of the error, 
where x a is again an approximate (e.g., computed) solution to the equation and 
|| || represents some vector norm. This error measure can be viewed as looking at 
the component by component size of the residual Ax n — b relative to \\b\\. This 
measure gives an accurate idea of the size of those components of the residual 
corresponding to the large components of b (i.e., those components of b which 
contribute significantly to ||t||). However, some small components of b may be 
poorly approximated by the corresponding components of Ax a without this error 
being reflected in the error measure \\Ax a - b\\/ ||fr||. Using \\b\\ as a comparison 
magnitude has masked the error in these small components. If one attempts to 
regain the lost information by looking at the size of each component of the residual 
vector relative to the magnitude of the corresponding component of b, one runs 
the risk of overestimating the error represented by components of the residual 
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corresponding to small components of b (to see the extreme case, imagine that 
some components of b are zero). 

This motivates the introduction of the notion of the data relative error. This 
is an error measure which is made on a component by component basis relative 
to comparison magnitudes which are based on all the data which went into the 
calculation of each component. In its ideal form, it is applicable to any arithmetic 
calculation producing a number which is theoretically supposed to be zero. For 
example, once an approximate solution x u to a linear system of equations Ax = b 
is known, the calculation Ax a — b of the residual has these properties. A calculation 
is made parallel to the calculation which produced the “should be zero” quantity. In 
the parallel calculation, each datum originally used is replaced by its absolute value, 
and all subtractions are replaced by additions. For the linear equation example, this 
parallel calculation is |-4||.r a | + |/>j, where all absolute values are taken component 
by component. 

Every individual scalar result x of the original calculation has its counterpart 
£ — 0 in the parallel computation, and (at least in exact arithmetic) |.r| < £. The 
magnitude of £ reflects the magnitude of the data which went into the calculation 
of x and gives an upper limit on what magnitude the answer could possibly have 
based on the form of the calculation. The calculation of £ does not take into 
account any cancellation, such as by the subtraction of nearly equal quantities, 
which might have occurred in the calculation of x. So, if £ > 0, the ratio |.r|/£ 
gives a reasonable idea of how close x is to 0 by comparison with the data from 
which it was calculated and what the calculation process could potentially do to 
that data. 

Unfortunately, the amount of programming necessary to calculate this ideal data 
relative error can be too large to be practical. For example, if .4, B, C, D'j 

and T [b) (^,A,B,C,D \ represent the results of computing the transfer function 

~ B + D in two different ways, then the quantity which is supposed 
to be close to zero is 

T [(l) (s, A, B, C, Dj - T {b) (,, I, B , C, Z >) . 

Computing the ideal data relative error would involve rewriting all of the software 
involved in calculations (a) and (b) to include parallel upper bound computations. 
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For purposes of Test 2 in this paper, a compromise was adopted. The compar- 
ison magnitude was taken from the computation 

c (si- a )' 1 \b\ + \d . 

The programming effort required to adapt the GP software to calculate (si — .4^ 
was moderate; the computational time to calculate the inverse was reasonable; 
and the results were reliable, even when the order of the matrix was 1445. The 
remainder of the calculation was routine. 

This comparison magnitude no longer has some of the properties of the ideal. 
For example, even if it were doubled to provide contributions for both transfer 
function calculations, the inequality |.r| < £ between an individual calculated 
number and its comparison magnitude no longer necessarily holds. This is because 

(si - Aj 1 is not an ideal comparison magnitude for (si - Aj . And 

i he compromise comparison magnitude does not necessarily refer to the actual 
algorithms used in computing either T< a) or T {b \ Despite this, the compromise 
comparison magnitude performed very well in assessing questions of accuracy in 
the present study. 


64 


References 


[1] Laub, Alan J.: Efficient Multivariable Frequency Response Computations. IEEE 
Transactions on Automatic Control , vol. AC-26, no. 2, 1981, pp. 407-408. 

[2] Laub, Alan J.: ALGORITHM 640. Efficient Calculation of Frequency Response 
Matrices from State Space Models. ACM Transactions on Mathematical 
Software , vol. 12, no. 1, 1986, pp. 26-33. 

[3] Maghami, Peiman G.; and Giesy, Daniel P.: Efficient Computation of Closed- 
Loop Frequency Response for Large-Order Flexible Systems. Journal of 
Guidance, Control, and Dynamics, vol. 19, no. 4, 1996, pp. 780-786. 

[4] Maghami, Peiman G.; Kenny, Sean P.; and Giesy, Daniel P.: PLATSIM: An 
Efficient Linear Simulation and Analysis Package for Large-Order Flexible 
Systems. NASA TP 3519, 1995. 

[5] Ogata, Katsuhiko: Modem Control Engineering. Prentice Hall, Second ed 
1990. 

[6] Golub, Gene H.; and Van Loan, Charles F.: Matrix Computations. The Johns 
Hopkins University Press, Second ed., 1989. 

[7] MATLAB High-Performance Numeric Computation and Visualization Software. 
Reference Guide. The MathWorks, Inc., Natick, MA, August 1992. 

[8] Lawson, C. L., Hanson, R. J.; Kincaid, D. R.; and Krogh, F. T.: Basic Linear 
Algebra Subprograms for Fortran Usage. ACM Trans. Math. Soft., vol. 5, no 3 
1979, pp. 308-323. 

[9] Lawson, C. L.; Hanson, R. J.; Kincaid, D. R.; and Krogh, F. T.: Algorithm 539. 
Basic Linear Algebra Subprograms for Fortran Usage [FI]. ACM Trans. Math. 
Soft., vol. 5, no. 3, 1979, pp. 324-325. 

[10] Dongara, J. J.; Moler, C. B.; Bunch, J. R.; and Stewart, G. W.: UNPACK Users’ 
Guide. SIAM, Philadelphia, PA, 1979. 

[11] Dongarra, Jack J.; Du Croz, Jeremy; Hammarling, Sven; and Hanson, 
Richard J.: An Extended Set of FORTRAN Basic Linear Algebra Subprograms. 
ACM Trans. Math. Soft., vol. 14, no. 1, 1988, pp. 1-17. 

[12] Dongarra, Jack J.; Du Croz, Jeremy; Hammarling, Sven; and Hanson, 
Richard J.: ALGORITHM 656. An Extended Set of Basic Linear Algebra Sub- 


65 


programs: Model Implementation and Test Programs. ACM Trans. Math. Soft., 
vol. 14, no. 1, 1988, pp. 18—32. 

[13] Dongarra, Jack J.; Du Croz, Jeremy; Hammarling, Sven; and Duff, Iain: A 
Set of Level 3 Basic Linear Algebra Subprograms. ACM Trans. Math. Soft., 
vol. 16, no. 1, 1990, pp. 1-17. 

[14] Dongarra, Jack J.; Du Croz, Jeremy; Hammarling, Sven; and Duff, Iain: 
ALGORITHM 679. A Set of Level 3 Basic Linear Algebra Subprograms: 
Model Implementation and Test Programs. ACM Trans. Math. Soft., vol. 16, 
no. 1, 1990, pp. 18-28. 

[15] Anderson, E.; Bai, Z.; Bischof, C.; Demmel, J.; Dongarra, J.; Du Croz, J.; 
Greenbaum, A.; Hammarling, S.; Mckenney, A.; Ostrouchov, S.; and Sorensen, 
D.: LAPACK Users’ Guide. SIAM, Philadelphia, PA, Second ed., 1995. 

[16] MATLAB Compiler User’s Guide. The MathWorks, Inc., Natick, MA, October 
1996. 

[17] Maghami, Peiman G.; Kenny, Sean P.; and Giesy, Daniel P.: PLATSIM: A 
Simulation and Analysis Package for Large-Order Flexible Systems ( Version 
2.0). NASA TM 4790, 1997. 

[18] Belvin, W. Keith; Maghami, Peiman; Tanner, Cheri; Kenny, Sean; and Cooley, 
Vic: Evaluation of CSI Enhancements for Jitter Reduction on the EOS AM-1 
Observatory. Dynamics and Control of Large Structures. Proceedings of the 
Ninth VPI&SU Symposium, L. Meirovitch, ed., Virginia Polytechnic Institute 
and State University, Blacksburg, VA, May 1993, pp. 255-266. 

[19] Giesy, Daniel P.; and Armstrong, Ernest S.: FREQ: A Computational Package 
for Multivariable System Loop-Shaping Procedures. NASA TM 4127, 1989. 

[20] Chiang, Richard Y.; and Safonov, Michael G: Robust Control Toolbox. For Use 
with MATLAB. The MathWorks, Inc., Natick, MA, August 1992. 

[21] Govaerts, W.; and Pryce, J. D.: Mixed block elimination for linear systems 
with wider boarders. IMA Journal of Numerical Analysis , vol. 13, no. 2, 1993, 

pp. 161-180. 

[22] Householder, Alston S.: The Theory of Matrices in Numerical Analysis. 
Blaisdell Publishing Company, New York, 1964. 


66 


Table 1. Time (seconds) to Calculate Frequency 
Response Function Using Original Implementation 


System 

size 0 

Freq. 

vector 

length 

Feed- 

forward 

present 

New 

MEX- 

code 

New 

M-code 

Old FOR- 
TRAN 
code 

Old 

M-code 

1/1 

24+39 

200 

No 

1.98 

11.18 

1.26 

7.57 

1/1 

24+39 

200 

Yes 

2.47 

13.50 

1.10 

7.62 

1/1 

24+39 

2000 

No 

24.78 

108.87 

10.08 

71.18 

1/1 

24+39 

2000 

Yes 

19.60 

135.20 

10.36 

71.92 

3/5 

184+39 

200 

No 

5.03 

18.53 

29.84 

287.93 

3/5 

184+39 

200 

Yes 

3.85 

14.18 

27.50 

290.25 

3/5 

184+39 

2000 

No 

46.08 

147.00 

213.69 

2818.78 

3/5 

184+39 

2000 

Yes 

35.77 

188.93 

200.09 

2830.32 

10/27 

1406+39 

200 

No 

60.40 

231.25 

5714.94 

49846.62 

10/27 

1406+39 

200 

Yes 

73.77 

227.40 

5756.99 

49199.20 

10/27 

1406+39 

2000 

No 

591.50 

1922.75 

24928.01 

- 

10/27 

1406+39 

2000 

Yes 

615.73 

1897.22 

24929.14 

- 


a In "m/n" in the first line, "m” is the number of inputs and "n” is the number of outputs. 

In "m+n" in the second line, "m" is the number of structural states and "n" is the number of controller states. 
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Table 2. Time (seconds) to Calculate Frequency 
Response Function Using Current Implementation 


Transfer 

function 

Number 
of system 
states 

New 

MEX- 

code 

New 

M-code 

Old FOR- 
TRAN 
code 

Old 

M-code 

GP code 

Tyr 

221 

14.26 

21.48 

37.64 

738.15 

142.52 

i V 

221 

21.74 

33.04 

50.58 

1232.54 

153.31 

Typr r 

221 

17.16 

27.89 

38.97 

764.57 

145.58 

Typr IC 

221 

30.00 

45.87 

53.93 

1278.97 

158.49 

Ter 

221 

14.47 

21.50 

35.88 

738.83 

142.71 

Tut 

221 

13.46 

20.23 

36.52 

731.16 

142.49 

T u( i 

221 

14.48 

22.03 

27.37 

368.50 

129.09 

Tyr 

1445 

27.27 

46.97 

5218.06 

- 

1413.02 

T yw 

1445 

83.86 

142.43 

5859.32 

- 

1550.30 

Typr r 

1445 

45.37 

87.68 

5436.88 

- 

1415.65 

TyPr w 

1445 

141.63 

257.97 

6052.90 

- 

1551.37 

T er 

1445 

27.34 

47.78 

5377.79 

- 

1407.49 

Tar 

1445 

22.05 

42.72 

5196.56 

- 

1422.89 

T \,d 

1445 

30.58 

51.21 

4664.98 

- 

1240.30 
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